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Successful operation of modern interconnected power 
systems depends largely on its ability to maintain stability 
during the transients. Transient stability evaluation assumes, 
therefore, greater significance in the planning, design and 
operation of large scale power systems. Power system stability 
is, in essence, a single concept; but for the purpose of ana- 
lysis, it is classified as 'transient stability' and 'steady 
state stability' . The former refers to the stability of the 
system under large disturbances such as faults on transmission 
network. Steady state stability refers to the ability of the 
system to remain in synchronism under small perturbations 
caused by random changes in the load. 


System security may be considered as the ability of a 
power system in normal state to undergo a disturbance without 



going into emergency condition. The static security assessment 
is widely used in modern energy control cehtres to predict the 
system security when subjected to a set of contingencies. It is 
assumed here that the system remains stable following a contin- 
gency. This assunption is relaxed in considering dynamic secu- 
rity. While the objectives and methodologies of dynamic security 
assessment are yet to be clearly defined, there is a general 
consensus that the system model is similar to that used in tran- 
sient stability studies just as static security assessment 
employs models used in load flov; studies. 

The transient stability analysis deals with nonlinear mathe- 
matical models describing the dynamic behaviour of the system. 
Digital simulation involving explicit evaluation of swing curves, 
is a versatile tool for stability assessment, but the method is 
time consuming and hence unsuitable for on-line applications. 
Direct methods based on transient energy calculations are recog- 
nised to be useful as screening tools in evaluating critical 
situations in the system. The transient energy margin (TEM) 
which is defined as the difference between the critical energy 
and, the energy at the clearing time has been suggested as an 
index for checking the stability of the system following a dis- 
turbance. This index also serves to indicate the severity of the 
disturbance. [ 



The energy-integral criterion, first proposed by Aylett, 
was later generalised in the use of Lyapunov functions by a 
number of researchers. The application of the Lyapunov method 
for transient stability analysis led to conservative results 
in the computation of the regions of stability. Later efforts 
using the energy function in the centre of inertia (COI) 
variables and the concept of controlling unstable equilibrium 
point (u.e.p) have removed much of the conservativeness of the 
results. Approximations in evaluating the uep's and the 
application of the Potential Energy Boundary Surface (PEBS) 
technique have greatly simplified the computation of stability 
regions . 

Hitherto, the application of direct methods using energy 
functions is based on reduced model of the system eliminating 
load buses. This is possible for constant inpedance type 
loads. Even if such a load representation is accurate, it 
will be desirable to preserve the structure of the network 
for on-line dynamic security assessment. The recent work on 
the structure preserving energy function (SPEF) (defined on 
structure preserving models) has removed this limitation and 
made it possible to consider arbitrary voltage dependent load 
models [l-3]. The problem of transfer conductances is also 
overcome by the use of structure preserving models. Some 
attempts have been made to include detailed machine models [4] 
where the results are conservative. 



While the concept of SPEF is undoutedly appealing, its 
application to on-line dynandc security assessment requires 
further developments. First of all, the computational require- 
ments with nonlinear load models are much more than with reduced 
models. This can be a severe limitation for on-line applica- 
tions. Secondly, there is a need to improve the system models 
used. Finally, the user must have the tools for the stability 
evaluation (both in planning and operation) that can be readily 
used for large systems. 

Another drawback of the present state of direct methods is 
the problem of accurate determination of the critical energy. 
While the PEBS method has overcome the computation of u.e.p, 
the accuracy of the method in multi-machine systems is governed 
by the fact that not all the kinetic energy (even with respect 
to centre of inertia) is responsible for system separation. The 
kinetic energy that is responsible for the system separation 
depends on the mode of instability and can be approximately cal- 
culated based on heuristic arguments [5]. However, an exact 
derivation of the energy function to be used in such cases is 
not available in the literature. 

This thesis is aimed at trying to overcome some of the 
limitations mentioned above. Specifically the objectives are : 

1) incorporation of the computation of SPEF in transient sta- 
bility program to simplify the stability evaluation. 



2) development of methods and efficient algorithms for the 
fast and accurate evaluation of transient stability of power 
systems thereby making the application of direct methods to 
on-line dynamic security analysis much more promising. 

3) development of SPEF for more realistic power system models 
for the evaluation of system stability. 

The major contributions of the thesis are listed below : 

i) An existing dynamic simulation program [6] is modified to 
include the computation of SPEF. The use of SPEF can speed 
up the calculation of critical clearing time (T._) or the 
stability limit thereby avoiding the time consuming simula- 
tion runs. Both classical and detailed machine models 
alongwith voltage dependent loads are used in the computa- 
tion of SPEF. 

ii) A new SPEF is derived considering the mode of instability 
where a single or a group of machines separate from the 
rest. This energy function differs from the previously 
derived energy function [3] in the kinetic energy term. 

The new energy function gives more accurate results. 

iii) Utilising the information on the mode of instability an 

algorithm is presented to speed up the confutation of TEM 
by solving algebraic equations. This procedure eliminates 
the need for determining the faulted trajectory and can 
make the on-line assessment of dynamic security using the 
energy function much more attractive. 



iv) A SPEF is derived to include the resistance of transmission 
lines which has not been done so far. Derivation is based 
on the assumption that conductance to susceptance ratios 
(G/B) of all the network elements are equal which permits 
network to be transformed into a lossless network [7] * This 
assumption is not found to be restricting as the numerical 
examples show that the use of an average value for g/B 
gives reasonably accurate results. 

v) A new SPEF is derived with realistic model of generators 
considering flux decay, transient saliency, AVR, exciter 
and damper winding on the quadrature axis. This is an 
extension of the work reported in [4]. However, the energy 
function derived here is constant on the post-fault 
trajectory with the time derivative as zero. This leads to 
more accurate results than that reported in [4] in the compu- 
tation of stability region using PEBS method. A simplified 
and general method for the calculation of SPEF is also 
presented where the use of classical model becomes a special 

C ct SO • 

The following is an outline of the work reported in this 
thesis : 

1) Chapter 1 introduces the various aspects of energy based 
direct methods for stability analysis of power systems and the 
application of these methods to on-line dynamic security assess- 
ment. The state of the art in this area is presented and the 
scope and objective of the thesis are outlined alongwith the 
chapterwise summary. 
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Chapter 2 deals with the modification of a dynamic simula- 
tion program by Podmore [6] to include the confutation of 
SPEF. The existing program provides time domain simulation 
of reduced povi/er system networks. It is written in FORTRAN 
and is composed of a ?AAIN routine and a number of. sub- 
routines beside two supporting programs for base case load 
flow and network reduction eliminating the load buses. The 
following modifications are incorporated in the existing 
program to develop the new transient stability program : 

A. The structure of the original network is retained to enable 
arbitrary voltage dependent load representation. This is 
done by 

i) replacing the supporting program for network reduction 
by the program YMAT to construct the network bus ad- 
mittance matrix. 

ii) modifying the subroutine NWSOL used for the network 
solution to include voltage dependent nonlinear loads. 

B. The program is augmented to compute SPEF. This is done by 
adding two new subroutines : 

i) PFSOL for solving the post-fault network equations at 
the end of each integration step to obtain voltages and 
their angles at all the buses to be used in the compu- 
tation of the energy function. Internal voltages of 
the generators to be used in the network solution are 
the same as that obtained by the integration of the 
machine equations. 



ii) TEF to evaluate the SPEF and its components for each time 

step of the network solution* 

The modified program can thus handle the power system sta- 
bility problem using i) time domain simulation and ii) direct 
assessment via the SPEF. The analysis of the program output 
during the faulted period can be used to predict the critical 
energy and the critical clearing time. After the initial deve- 
lopment, the transient stability program was tested with three 
realistic power system examples. The validation test confirms 
that direct methods provide results comparable to that of simula- 
ti on. 

3) Chapter 3 is devoted to the systematic development of a 
mathematical model for the direct stability evaluation of multi- 
machine power systems based on the assumptions that i) instabi- 
lity results from a single machine or a group of machines separa- 
ting from the rest and ii) the generators in each group swing 
together. A new energy function is then derived and it is shown 
that the function has a time derivative which is zero. Structure 
preserving model of the system is used and the generators are 
represented by the classical model. Loads are modelled as 
arbitrary functions of respective bus voltage. The system equa- 
tions and the energy function are formulated in COI variables. 
Critical energy is evaluated using the PEB3 method. The proce- 
dure has been applied to three system examples considered in 
Chapter 2. The results indicate that i) the inclusion of the 



mode of instability in the energy function gives consistently 
a better estimate of the ci^itical clearing time conpared to 
that obtained using the energy function of ref. [3] which does 
not consider the mode of instability and ii) the effect of 
voltage dependent active power loads on the region of stability 
is quite significant and thus cannot be ignored. The voltage 
dependent active power loads introduce path dependent term in 
the energy function which is computed by numerical integration 
using trapezoidal rule. 

4) An efficient method for assessing the dynamic security of 
power systems on-line based on the concept of TEM» is presented 
in Chapter 4. The system model and the energy function of 
Chapter 3 are again used for the development of the method. A 
deterministic estimation of ^^9 transient energy margin invol- 
ves the computation of the energy function i) at criticality 
and ii) at clearing. 

It is shown that for a given mode of instability (as des- 
cribed in Chapter 3) f the potential energy reduces to a fun- 
ction of single independent variable corresponding to the rotor 
angle(s) of the machine(s) going out-of-step with respect to 
the remaining machines* It is also shown that the critical and 
the clearing energies can be obtained by the solution of alge- 
braic equations thereby completely eliminating the simulation of 
the faulted system. This leads to a substantial reduction in 
the computation required for stability evaluation using SPEF. 



The proposed method is illustrated by considering the 
three system exanples given in Chapters 2 and 3. The predicted 
value of the critical energy and that at clearing are in close 
agreement with those obtained using the faulted trajectory simu- 
lation. The proposed method can also be used for accurate rank- 
ing of the contingencies as illustrated in 10-machine system 
example. 

5) Structure preserving energy functions have so far been 
developed without accounting for the transmission line resis- 
tances. An attempt is made in Chapter 5 to develop a valid 
energy function including the line resistances. The method is 
based on a network property which states that a network having 
the same G/b ratio for all its elements is equivalent to a loss- 
less network with a new set of power injections [7]. The use of 
the proposed energy function is illustrated with numerical exam- 
ples. Results indicate that i) inclusion of line resistances 
improves transient stability in terms of a larger region of 
stability and a higher value of critical dclearing time and 
ii) predicted critical clearing time obtained using the average 
value of G/b for the network elements conpares well with that 
obtained by digital simulation using the actual resistance of 
these elements. The component of the energy function accounting 
for the transmission line resistances is path dependent and is 
evaluated by trapezoidal rule. 



6) The principal aim of Chapter 6 is to extend the work 
reported in ref. [4] considering more detailed models of the 
generators. The topological energy function of ref. [4] has a 
time derivative which is negative semi-definite. Considering 
the same model as in [4], a new SPEF is developed, the time 
derivative of which is zero. Moreover, this model is improved 
upon by the addition of a damper winding on the quadrature axis 
and a new energy function is derived. The time derivative of 
the energy function is again shown to be zero. Simpler expre- 
ssions for the energy functions are then derived thereby 
increasing the speed of calculation. Investigations are carried 
out on a 4-machine, and a- lO-machine' system examples. Critical 
clearing times obtained by direct method and by simulation are 
compared and are found to be in close agreement. The results 
show that i) the variation of the field and the generator vol- 
tages have predominant effects on the system energy during a 
transient and ii) the inclusion of transient saliency, AVR and 
damper winding is beneficial in inproving the transient 
stability. 

7) A brief review of the major contributions of this thesis 
and the suggestion for further work form the seventh and the 
concluding chapter. 
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CHAPTER 1 


INTRODUCTION 


1.1 GENERAL 

There has been a phenomenal growth, in size and complexity, 
of electric power systems to match the increasing demand for 
electric energy. Successful operation of these modern inter- 
connected systems depends largely on their ability to maintain 
stability during transients. Transient stability analysis, 
therefore, assumes greater significance in the planning, design 
and operation of large scale power systems. Power system sta- 
bility may be defined as that property of the system which 
enables it to respond to a perturbation from a normal operating 
point so as to return to a condition where its operation is 
again normal [l-8]. The main criterion for stability is that 
the system -maintains synchronism at the end of the transient 
period. 

Power system stability is, in essence, a single concept; 
but for the purpose of analysis, it is classified as 'transient 
stability’ and ' steady state stability' depending upon the 
nature and the magnitude of the disturbance. Transient stabi- 
lity studies are aimed at determining system stability follow- 
ing major disturbances. These disturbances may be due to 
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transmission system faults* sudden injection or rejection of 
a large load, tripping of a large generating unit or line 
switching. Three phase short circuits followed by line 
clearing are the most severe of all these disturbances. Usually* 
disturbances alter the system at least temporarily so that sub- 
sequent steady state operation will be different from that 
prior to the disturbance* There is then the necessity for that 
altered system to be stable in its new steady state. 

Steady state stability refers to the ability of the system 
to remain in synchronism under small perturbation caused by 
random changes in load. Decrease in transmission system capa-* 
bilities or inadequately designed regulating devices (such as 
AVRs) can cause steady state instability. 

System security may be considered as the ability of a 
power system in normal state to undergo a disturbance without 
going into emergency condition. The static security assessment 
is widely used in modern energy control centres to predict the 
system security when subjected to a set of contingencies. It is 
assumed here that the system remains stable following a contin- 
gency. This assumption is relaxed in considering dynamic secu- 
rity. While the objectives and methodologies of dynamic : 
security assessment are yet to be clearly defined, there is a 
general consensus that the system model is similar to that used 
in transient stability studies just as static security assess- 
ment employs models used in load • flow .studies * 
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Methods of Analysis 

The transient stability analysis of power systems deals 
with nonlinear mathematical models describing the dynamic 
behaviour of the systems. These models consist of a set of 
differential equations representing the dynamics of generating 
units (synchronous machines and their controls) and a set of 
algebraic equations describing the interconnections through the 
transmission network. Various approaches to the stability pro- 
blem are : 

i) Digital simulation 

ii) Direct methods 

iii) Pattern recognition 

Digital simulation techniques involve the numerical 
integration of the system differential equations to obtain 
explicit time domain solutions. The numerical integration 
methods used for stability analysis, although very effective 
in handling sophisticated models, are time consuming and hence 
unsuitable for on-line applications. 

To overcome this limitation, energy-based direct methods 
for stability analysis have been developed. Computational 
effort is significantly reduced since only the faulted system 
is integrated for a short time till the clearing of the fault. 
As a result of continuous research and the progress made, these 
methods are now recognised to be useful as screening tools in 
assessing critical situations in the system [9]. 
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Pattern recognition [10-12] is another approach aimed at 
overcoming the high computational requirements of on-line 
transient stability studies. A large number of off-line sta- 
bility studies are performed to form a ’ training set' and 
certain important features are selected. A classifier compares 
the actual operating conditions with the training set and, on 
the basis of this comparison, classifies the existing state as 
either secure or insecure. 

Direct methods for stability evaluation involve the com- 
parison of the energy attained by the system during the fault- 
on period with the critical energy. The transient energy 
margin (TEM) which is defined as the difference between the 
critical energy and the energy at the clearing time, has been 
suggested as an index for checking the stability of the system 
following a disturbance. This index also serves to indicate j 

the severity of the disturbance. Normalised transient energy I 

j 

margin which may be defined as the excess potential energy j 
capability (difference between the critical energy and the 

! 

potential energy at fault clearing) of the system, expressed I 

as the ratio of the kinetic energy at the instant the fault is i 

i 

cleared, also serves as an alternative measure for inferring | 

i 

system stability or instability. ! 

The application of direct methods to transient stability 
evaluation and on-line dynamic security assessment of power 
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systems using structure preserving energy functions (SPEF) is 
investigated in this thesis. 

1.2 DIRECT METHODS 

As stated earlier, transient stability analysis is 
concerned with the system’ s ability to remain in synchronism 
following major disturbance such as a three phase fault at the 
terminals of generators. During transients the dynamics of the 
system can be described by the following set of differential- 
algebraic equations which are functions of the system variables 

i - I. (2£» l) (l*la) 

(1.1b) 

i 

The functions f and g will be different for the faulted and the; 
post-fault networks because of change in network configuration.: 

Transient stability study using numerical integration 

i 

methods involves explicit time domain solutions of Eqs. (1.1). ; 

1 

These methods are, therefore, expensive in terms of computation! 
requirement. For this reason, search for direct methods, that ; 

i 

is, methods to determine stability independently of the time ; 
solution of Eqs. (1*1), has continued. The classical equal 
area criterion is such a direct method for single- or two- 
machine systems. 



The post-fault system has an initial state decided by 
the state at the end of the fault period. Also, the post- 
fault system has a stable equilibrium point (sep) which can 
be obtained by solving the set of equations in (1.1) corres- 
ponding to the post-fault network. The fundamental problem 
of transient stability is to decide whether the post-fault 
system, starting in the initial state, will coverge to the sep 

There is a region of initial conditions in the state 
space from which trajectories converge to the sep. This 

region is called the region of stability. If the initial post 

fault state lies within this region, then the system is stable 

Direct methods attempt to compute this region through the com- 

putation of a Lyapunov or energy function which is a positive 
definite function of the system variables as shown below : 

W = W (x, j) (1.2) 

and whose time derivative along the post-fault trajectory is 
negative (semi-) definite. The region of stability surround- 
ing the sep can be defined by 

W (x, y) < (1»3) 

where W is the critical value of W. In testing transient 
cr 

stability following a fault, the value of W is monitored at 
the time of clearing (t^j^) and compared with If 
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then the system is considered stable. The evaluation of WCt^j^) 
requires the solution of Eqs. (1.1) corresponding to the faulte 
system configuration. This is usually obtained by simulation. 

The critical clearing time (t ) which is of interest in 

C IT 

stability studies, is defined by 


W(t 




(1.5) 


1.3 LITERATURE SURVEY 
1.3.1 General 

The use of equal area criterion is the first application 
of Lyapunov's direct method to power system transient analysis. 
However, the systematic application of this method to the 
stability problem of realistic power systems is due to Gless 
[13] and El-Abiad and Nagappan [14] . Since then there has been 
growing interest in the application of direct methods to tran- 
sient stability analysis. A summary of the research efforts is 
well documented in the survey papers by Ribbens-Pavella [15,16] 
Fouad [17], Pai [18] , Varaiya [19] and in a recent monograph on 
the subject by Pai [l]* These investigations have essentially 
been directed along two distinct, although related lines : 

i) the construction of suitable Lyapunov functions with 

improved system models to allow for flux decay and regula- 
ting devices such as AVR, governor, etc. [20-26], 
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ii) development of fast and accurate algorithms for implemen- 
tation of the method to on-line transient stability stu- 
dies [1, 27-33]. 

1.3.2 Lyapunov or Energy Functions 

The systematic construction of Lyapunov functions was 
attempted by Willems [20, 21] and Pai et al [22]. Detailed 
models of machines were considered by Pai and Rai [23] for 
stability analysis of synchronous machines using Lyapunov 
Popov criterion. Results obtained with models higher than 
the second order representation of machines appear to be 
conservative in addition to the difficulties involved in the 
construction of an effective Lyapunov function. This motiva- 
ted Sasaki [24] to approximately include the field flux decay 
in the transient stability analysis. A Lyapunov function 
based on Moore-Anderson Theorem considering flux decay in 
generators was constructed in ref. [25]. In their recent 
work, Kakimoto et al [26] have incorporated more realistic 
model for generators including AVR and PSS in investigating 
transient stability by direct methods. However, the Lyapunov 
function in [25] does not include the effect of AVR and gives 
inconsistent results for the predicted regions of stability. 

Energy functions are possible Lyapunov functions for 
conservative systems. The various components of the energy 
function can be given a physical interpretation. Energy based 
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concepts for power system transient stability analysis was 
first proposed by Magnusson [34] and then by Aylett [35] . The 
energy-integral criterion was later generalised in the use of 
Lyapunov functions. However, for about last one decade, there 
has been renewed interest in the use of energy function for 
stability analysis, particularly dynamic security assessment. 

The centre of inertia (COI) is an important variable in 
power system analysis. It is defined as the inertia weighted 
average of all the rotor angles [36, 37]. Reference [36] shows 
that for transient stability evaluations, the kinetic energy of 
the system must be calculated by utilising the relative veloci- 
ties with respect to CX)I . In ref. [30], COI variables are used 
to formulate the system equations and to develop an energy fun- 
ction with reduced systems. 

In their recent work, Pai et al [38, 39] have exploited a 
physically based decomposition technique to perform direct sta- 
bility analysis of power systems using energy functions. For- 
mulating the system equations in a singularly perturbed form, 
they use the slow coherency based grouping technique [40,41] to 
decompose the system into r areas. This grouping is then used 
to develop i) slow energies associated with the centres of 
inertia of the areas, ii) fast energies associated with inter- 
generator oscillations within an area and iii) slow-fast ener- 
gies associated with both the slow and the fast variables. It 
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has been shown that depending on the fault location, only a. 
few areas are disturbed and it is possible to approximate the 
system energy function by the fast energies of the disturbed 
areas, global slow energy and the slow-fast energy. 

VVinkelman et al [42] have proposed a Hamiltonian based 
approach to study the transient stability problem including 
flux decay and voltage regulation. They have developed a 
generalised potential energy function including the effects 
of non-conservative forces acting on the system during tran- 
sients. The total energy of the post-fault system has been 
shown to remain constant. 

Transfer Conductances 

Hitherto, the application of direct methods using energy 
functions is based on reduced model of the system eliminating 
the load buses. This is possible for constant impedance type 
loads. Such a load representation introduces transfer con- 
ductances in the nodal admittance matrix (corresponding to the 
internal nodes of the generators) . Aylett [35] has shown that 
the presence of transfer conductances makes some of the terms 
in the Lyapunov function path dependent. He, therefore, 
suggested the omission of these terms or for the sake of 
consistency, the omission of all the transfer conductance terms. 
Later, Ribbens-Pavella [43] successfully developed two valid 
Lyapunov functions following the procedure of Aylett. 
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Uemura et al [44] handled the problem of transfer 
conductances by making approximations in the integrand which 
was once considered as a constant and another time as a linear 
function. Pai et al [45] succeeded in constructing valid 
Lyapunov functions with transfer conductances for a tv;o-machine 
system. They claimed that the generalisation of their approach 
to multi -machine systems was straightforward. This prompted 
Gudaru [46] to derive a Lure type Lyapunov function for the 
multi -machine systems with transfer conductances following the 
Anderson procedure. Henner [47] and Willems [48] have, however, 
shown. that the generalised Lyapunov function derived by him is 
not valid since he fails to satisfy the positive real matrix 
criterion . 

Later, Pai and Varwandkar [49] used the same non-lineari- 
tios as Gudaru and the Lyapunov function of Pai and Murthy [45] 
neglecting the transfer conductances and succeeded in generating 
the Lyapunov function by introducing approximations in the model 
for the multi -machine case. 

Kakimoto et al [26] considered the effect of transfer con- 
ductances by adding numerical correction to the Lure type 
Lyapunov function under the assumption that the system traje- 
ctories are not so much different in the relative angular space 
between a sustained fault and the critically cleared fault and 
obtained two modified Lyapunov functions. 
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Athay et al [30] computed the critical value of the energy 
function using a linear path of integration while for evaluating 
the energy function along the faulted trajectory, they used the 
actual path to integrate the non-integrable terms via the tra- 
pezoidal rule. 

Guindi and Mansour [50] approximated the relative velocity 
by a constant factor and succeeded in constructing a valid 
Lyapunov function. 

Structure Preserving Energy Functions (SPEF) 

Another disadvantage of using the reduced system model is 
that it masks the topology of the system network to be investi- 
gated. Beside, it does not allow one to understand and inter- 
pret the, energy shifts among system components during a distur- 
bance. The problem of transfer conductances due to loads can 
be avoided by an approach suggested by Pai et al [51], using 
distribution factors to reflect the effect of loads as current 
injections at the generator internal buses. Later efforts using 
structure preserving energy functions (SPEF) have greatly helped 
in removing this limitation and made it possible to consider 
more realistic load models. 

The concept of structure preserving energy functions to be 
used for transient stability evaluation is due to Bergen and 
Hill [52] who first proposed a system model retaining the iden- 
tity of the load buses. They model the load buses as P-V type 
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buses and allow the active powers of load to be frequency 
dependent; but their ’ topological’ energy function does not 
depend on frequency dependence. In addition, the assumption of 
constant bus voltages is not valid, since these voltages flu- 
ctuate considerably during a fault. Their approach has been 
improved upon by Narasinghamurthy and Musavi [53] by consider- 
ing voltage dependent reactive power loads in developing the 
SPEF. However, the real power at each load bus is assumed to be 
constant. An improved load model has been proposed by Padiyar 
and Sastry [54] and Ribbens-Pavella et al [55] by considering 
the real and the reactive power of load to be arbitrary fun- 
ctions of respective bus voltage. In their recent work, 
Ribbens-Pavella et al [56] have presented an unified approach 
for modelling the loads by allowing the real power load at each 
bus to be both voltage and frequency dependent and the reactive 
power load to be voltage dependent. Attempts have also been 
made to include detailed machine models in the development of 
SPEF. Tsolas et al [57] developed the structure preserving 
energy function considering flux decay and transient saliency 
in machines using constant active and reactive power loads. In 
ref. [58], Padiyar and Sastry developed a ’topological* energy 
function considering flux decay, transient saliency, exciter and 
AVR allowing the loads to be voltage dependent. However, their 
energy function has a time derivative which is negative semi- 
definite and the results obtained appear to be conservative. 
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1.3.3 Computation of Stability Regions 

Direct methods for power system transient stability 
evaluation involve an accurate determination of the region of 
stability or the critical energy. In the past, the minimum 
value of the Lyapunov function evaluated at all possible unsta- 
ble equilibrium points (uep) was chosen [l]. This procedure 
yields very conservative results for the regions of stability. 
Besides, the computational burden increases prohibitively for 
large systems. This prompted Prabhakara et al [27] to suggest 
simplifications in the determination of these points. Ribbens- 
Pavella [31] has shown with various numerical examples that 
consideration of the common modes of instability reduces the 
number of uep' s to be evaluated to 2m, m is the number of 
machines. In ref. [30], Athay et al have shown that consi- 
deration of the controlling or the relevant uep substantially 
eliminates the limitations of the method due to conservative- 
noss of results. 

The Potential Energy Boundary Surface (PEBS) method has 
overcome the need for the computation of uep' s [29, 30]. PEBS 
is defined as the surface formed in angle space by the points 
corresponding to the first maxima of the potential energy with 
respect to the stable equilibrium point, where its directional 
derivative is zero. The uep's (where the gradient of potential 
energy is zero) lie on the PEBS. The critical energy is deter- 
mined as the potential energy at the instant when the faulted 



X5 


trajectory crosses the PE3S . It is assumed that the critically 
cleared trajectory passes very close to the controlling uep. 

The results obtained using the PEBS method are, usually, in 
good ag|reement with those obtained by the controlling uep 
method . 

Although the PEBS method succeeds in avoiding the need to 
compute the uep’s, the accuracy of the method in multi -machine 
system is governed by the fact that not all the kinetic energy 
is responsible for system separation. The kinetic energy that 
is responsible for system separation depends on the mode of 
instability and can be approximately calculated based on heuri- 
stic arguments [59]* The concept of generalising the equal 
area criterion for multi -machine systems, is proposed in ref. 
[32] by Fouad and Vittal. This is used to determine the con- 
trolling uep and the critical energy for the disturbance under 
investigation [33]. Michel et al [60] have used the transient 
energy of individual machines to study the stability problem. 
They claim that the use of individual machine energy functions 
rather than the overall system energy helps in correctly pre- 
dicting the mode of instability and the critical clearing time. 

1.3.4 Application to on-line Dynamic Security Assessment 

The problem of dynamic security assessment in an on-line 
environment remains largely unsolved and there is not much 
literature available. In refs. [32, 61], Fouad et al attempt 
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to assess the transient stability of power systems through the 
application of transient energy margin (TEM) . With the help 
of numerical examples, they demonstrate that the use of TEM 
provides information on the weak links in the system and helps 
in ranking three phase and single line- to-ground faults. 
Ribbens-Pavella et al [^2] derive a transient stability index 
(TSI) based on the maximum initial relative acceleration of the 
system. They then use this index to assess the robustness of 
a power system subject to a disturbance, without performing any 
disturbance simulation. However, their approach is heuristic 
in nature. Pai et al [63, 64] propose the use of energy mar- 
gins for economic load dispatch with dynamic security conr .. 
straints. They exploit the sensitivity of the energy margins 
to determine maximum load capability or simultaneous inter- 
change capability. In their recent work, Fouad et al [65] 
suggest a technique incorporating the sensitivities of the un- 
normalised transient energy margin to speed up the transient 
energy function method in on-line dynamic security assessment 
of power s^ystems. 

There have been no attempts for the application of SPEF 
for on-line dynamic security assessment, reported in the 
literature. 
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1.4 OBJECTIVES mU SCOPE OF THE THESIS 

While the concept of structure preserving' energy function 
is undoubtedly appealing, its on-line application to dynamic 
security assessment requires further developments. First of 
all, the computational requirements with nonlinear load models 
are much more than with reduced models. This can be a severe 
limitation for on-line applications. Secondly, there is a need 
to improve the system models used. Finally, the user must have 
the tools for stability evaluation (both in planning and opera- 
tion) that can be readily used for large systems. 

This thesis is aimed at trying to overcome some of the 
limitations mentioned above. Specifically, the objectives are; 

i) incorporation of the computation of SPEF in transient 
stability program to simplify the stability evaluation. 

ii) development of methods and efficient algorithms for the 
fast and accurate evaluation of transient stability of 
power systems thereby making the application of direct 
methods to on-line dynamic security analysis much more 
promising . 

iii) development of SPEF for more realistic power system 
models for evaluation of system stability. 
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1.5 OUTLINE OF THE THESIS 

The chapterwise summary of the work done in this thesis 
is given below ; 

1) Chapter 2 deals with the modification of a dynamic simula- 
tion program by Podmore [66] to include the computation of SPSF» 
The existing program provides time domain simulation of reduced 
power system networks. It is written in FORTRAN and is com- 
posed of a MAIN routine and a number of subroutines beside two 
supporting programs for base case load flow and network redu- 
ction eliminating the load buses. The following modifications 
are incorporated in the existing program to develop the new 
transient stability program : 

A. The structure of the original network is retained to enable 
arbitrary voltage dependent load representation. This is 
done by 

i) replacing the supporting program for network reduction 
by the program YMAT to construct the network bus admi- 
ttance matrix 

ii) modifying the subroutine NWSOL used for network solu- 
tion to include voltage dependent nonlinear loads. 

B The program is augmented to compute SPEF. This is done by 
adding two new subroutines : 
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i) PFSOL for solving the post-fault network equations at 
the end of each integration step to obtain voltages 
and their angles at all the buses to be used in the 
computation of the energy function* Internal voltages 
of the generators* to be used in the network solution 
are the same as that obtained by the integration of 
the faulted system equations. 

ii) TEF to evaluate SPEF and its components for each time 
step of the network solution. 

The modified program can thus handle the study of the 
power system transient stability problem using a) time domain 
simulation and b) direct assessment via the SPEF* The analysis 
of the program output during the faulted period can be used to 

\ 

predict the critical energy and the critical clearing time* 

After the initial development, the transient stability program 
is tested with three realistic power system examples* The 
validation test confirms that direct methods provide results 
comparable to that of simulation* 

2) Chapter 3 is devoted to the systematic development of a 
mathematical model for the direct stability evaluation of multi* 
machine power systems based on the assumptions that i) instabi- 
lity results from a single machine or a group of machines sepa- 
rating from the rest and ii) the generating units in each group 
swing together* A new energy function is then derived and it 
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is shown that the function has a time derivative which is 
zero along the post-fault trajectory. The generators are 
represented by the classical model. The machine damping and 
transmission line losses are ignored. Loads are modelled as 
arbitrary functions of respective bus voltage. The system 
equations and the energy function are formulated in COI varia- 
bles. Critical energy is evaluated using the PEBS method. 

3) An efficient method for assessing the dynamic security of 
power systems on-line, based on the concept of TEM, is presented 
in Chapter 4. The system model and the energy function of 
Chapter 3 are again used for the development of the method. A 
deterministic estimation of the transient energy margin invol- 
ves the computation of the energy function i) at criticality 
and ii) at clearing. 

It is shown that for a given mode of instability (as des- 
cribed in Chapter 3) , the potential energy reduces to a fun- 
ction of single independent variable corresponding to the rotor 
angle of the machine going out-of-step with respect to the 
remaining machines. It is also shown that the critical. and 
the clearing energies can be obtained by the solution of al- 
gebraic equations thereby completely eliminating the simulation 
of the faulted system. This leads to a substantial reduction 
in the computation required for stability evaluation using 


SPEF 
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The proposed method is illustrated by considering the 
three system examples given in Chapters 2 and 3. 

4) Transmission lines have so far been considered to be 
lossless in the derivation of structure preserving energy fun- 
ctions. An attempt is made in Chapter 5 to include the line 
resistances in the development of a valid energy function. The 
method is based on the network property that a network, all 
the elements of which have the same G/b ratio, is equivalent 
to a lossless network with a new set of power injections [67], 
An energy function is derived based on the transformed network. 
It is shown that the time derivative of the function on the 
post-fault trajectory vanishes. The use of the proposed energy 
function is illustrated with numerical examples. 

5) The principal aim in Chapter 6 is to develop a structure 
preserving energy function including detailed generator models. 
This work is based on an earlier work by Padiyar et al [58 ] 
wherein they have derived a ’ topological' energy function 
incorporating the effect of flux decay, transient saliency, 
exciter and AVR. The function in [58] has a time derivative 
which is negative semi -definite along the post-fault trajectory 
and gives pessimistic estimates of the region of stability. An 
attempt has been made in this chapter to remove these limita- 
tions. The generator model used in [58] is augmented by the 
addition of a damper winding on the quadrature axis. A new 
SPEF corresponding to the improved model is then derived. The 
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time derivative of the energy function is shown to be zero 
along the post-fault trajectory. A simpler and yet general 
expression for the energy function is then derived where the 
use of the classical model becomes a special case. The new 
energy function is tested on a 4-machine and a 10-machine system 
examples. 

6) The summary of the research results in this thesis and 
suggestions for further work form the seventh and the conclu- 
ding chapter. 



CHAPTER 2 


POWER SYSTEM DYNAMIC SIMULATION INCORPORATING 

COMPUTATION OF STRUCTURE PRESERVING ENERGY 

FUNCTIOJ 


2.1 INTRODUCTIOl 

Transient stability evaluation using digital simulation 
is well developed. There are efficient transient stability 
programs which incorporate detailed generator models and 
voltage dependent load models. These programs are usually 
used for planning studies. 

The use of direct methods can permit on-line evaluation 
of transient stability for security analysis. Even for off- 
line planning studies » direct methods can help in reducing 
the time of computation by reducing the number of cases that 
need to be studied in detail, in other words, acting as 
screening tools [9]. This is applicable even though sim- 
plistic models (classical models) of generators are used in 
direct stability analysis. 

A major factor that has prevented the extensive use of 
direct methods by system planners is the nonavailability of 
effective computer programs for the computation of energy 
functions. Actually, modifying the existing programs can 
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result in a transient stability program which can help a 
user to apply direct methods or simulation for studying a 
particular case. The importance of this aspect has been 
appreciated by researchers [68]. It is also possible to 
use energy functions to provide fast information on the 
system stability characteristics, particularly in large 
systems where the handling of voluminous results, even in 
the form of swing curves, can be cumbersome. 

Structure preserving energy function (SPEF) has the 
advantage of handling voltage dependent loads and even 
detailed generator models [57,58]. The objective of this 
chapter is to incorporate the computation of SPEF in an 
existing stability program and report on the results obtained 
using the augmented program. Power system dynamic simulation 
program by Podmore [66] is available and is well suited for 
research applications. This program enables detailed repre- 
sentation of generators with their controllers. The only 
restriction is the representation of load as constant impe- 
dances which permits network reduction. This limitation of 
the program is removed by rewriting some of the subroutines 
used in the program. The revised program can handle any 
arbitrary voltage dependent load models. 

The augmented program gives additional output in terms of 
the energy function and its various components. By simulating 
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the power system only during the faulted period, it is 
possible to predict the critical energy using the potential 
energy boundary surface (PEBS) method, from the analysis of 
energy function. The application of the augmented program for 
stability evaluation is presented using case studies on three 
different systems (7-machine, 10-mac hine and 13-machine systems) 
These results are also validated by simulations using the same 
program. 

2.2 POWER SYSTEM DYNAMIC SIMULATION 
2.2.1 General 

The equations representing the dynamic behaviour of a 
power system can be conveniently divided into two groups : 

i) Equations describing the dynamic behaviour of the machines 
and their controllers. These consist of a set of differential 
equations of the form : 

x=f(x,y) (2.1) 

ii) Equations describing the steady state behaviour of the 
network and generator armature circuits. These consist of 
algebraic equations of the form : 

g(x,y) = 0 (2.2) 

Variables x are defined as the system state variables. 
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Variables y are referred to as auxiliary variables and are 
associated with the network. Generator speed and angular 
position are examples of state variables, whereas bus voltages 
and bus angles are examples of auxiliary variables. 

The structure of the network equations (2.2) may change 
in time due to network changes such as fault inception, fault 
clearing, line switching etc. 

Any dynamic simulation program has to solve Eqs. (2.1) and 

(2.2) . Numerical integration is generally used to solve the 
nonlinear differential equations (2.1). The network equations 

(2.2) reduces to a linear equation if the loads are assumed 
to be constant impedances. In this case, y can be solved 
explicitly in terms of x and substituted in Eq. (2.1). 

2.2.2 Details of Existing Simulation Program [66] 

This program is written in FORTRAN and is composed of a 
IvlAIN routine and a number of subroutines. The program requires 
the support of two other programs : 

i) A power system load flow program 

ii) A program for network reduction eliminating the load buses. 

In the simulation program, loads are assumed to be constant 
impedances, which enables the load buses to be eliminated using 
network reduction. Only the generator buses are retained and 
one generator is assumed to be connected to each bus. 
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The simulation program is structured as follows : 

i) MAIN routine - for data initialization, input and 
providing overall control of the program. 

ii) NWSOL - 9 subroutine for solving the armature and the 
network equations. This subroutine calculates the 
terminal voltages and the bus current injections 
(generator outputs) of the network. 

iii) Equipment Subroutines ~ a library of subroutines for 
modelling various types of excitation systems and 
turbine-governor systems. Each equipment subroutine 
consists of two parts, one part calculates the initial 
values of the state variables and is executed once at 
the start of the program, the other calculates the time 
derivatives of the state variables for each integration 
time step. 

The various synchronous generator representations are 
modelled with a single subroutine GENl. The routine receives 
the terminal voltage and current as input and provides the 
initial conditions of the generator state variables. The 
synchronous machine model includes flux decay, AYR, transient 
saliency and a damper winding on the quadrature axis. By 
choosing the parameters appropriately, it is possible to 
represent simpler models of the machine (say, classical model). 
The machine model is presented in APPENDIX A. 
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The excitation system is modelled using three subroutines! 

AVRl : IEEE Type 1 representation of rotating excitation 
system 

AVR2 : IEEE Type is representation of static excitation system 
AVR3 : Rotating excitation system with an auxiliary stabiliser. 

The turbine governor modelling is done by the following 
subroutines : 

TURl ; Three stage steam turbine (HP, IP and LP) and governor 
TUR2 : Hydraulic turbine and governor. 

iv) INT - an integration subroutine which calculates the 
state variables for a new time interval as a function of 
the state variables and state variable derivatives for 
the present and the preceeding intervals. The integratin'' 
routine is based on the formula : 

^t+ it = it + (it - it- At> -St (2.3) 

where x^ and are the state variables and their time 
derivatives at time indicated by the subscript. At is 
the time step. 

v) OUTPUT and PLOT - these subroutines respectively give 
printouts and plots of the variables calculated during 
the simulation. 
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During the calculation it is necessary to account for the 
interconnections between the generator, excitation system and 
turbine governor models. Fig. 2.1 shows the relationship 
between the various subroutines for the calculation of the 
initial conditions. In this case, the field voltage and the 
mechanical power are calculated by the generator subroutine and 
consequently, this is executed before the excitation system 
and turbine-governor subroutines. The relationship which 
exists between the subroutine input and output variables during 
the step by step simulation is shown in Fig. 2.2. In this case, 
the new values of field voltage andl mechanical power are res- 
pectively calculated by the excitation system and turbine- 
governor subroutines. Consequently, the generator subroutine 
is now executed after the excitation system and turbine- 
governor subroutines. 

2.3 modification OF THE DYNAiVlIC SIMULATION PROGRAM 
2.3.1 General 

The following modifications are incorporated in the 
existing’ program to develop the new transient stability 
program : 

Incl usion of Voltage Dependent Loads 

Topology of the original network is preserved to enable 
loads to be modelled as arbitrary functions of respective bus 




FI6.2.1 RELATIONSHIP OF SUBROUTINE VARIABLES FOR 
CALCULATING INITIAL CONDITIONS. 



FIG. 2.2 RELATIONSHIP OF SUBROUTINE VARIABLES DURING 
STEP-BY-STEP SIMULATION. 
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voltage* This is done by 

i) replacing the supporting program for network reduction by 
the program VMAT to construct the network bus admittance 
matrix 

ii) modifying the subroutine NWSOL used for the network solu- 
tion to include voltage dependent nonlinear loads. 

An algorithm given in rof. [69] is used for solving the 
combined network and generator armature equations. This 
handles both saliency and nonlinear loads. The method requires 
iteration at each integration time step. The details of the 
algorithm are given in APPENDIX B. 

Inclusion of SPEF 

The program is augmented to include the confutation of 
SPEF. This is done by adding two new subroutines : 

i) PFSOL for solving the post-fault network equations at the 
end of each integration step to obtain voltages and their an- 
gles at all the buses to be used in the computation of the 
energy function. Internal voltages of the generators to be 
used in the network solution are the same as that obtained by 
the integration of the machine equations. This routine is to 
be used only during the faulted period. 

ii) TEF to evaluate the SPEF and its components for each time 
step of the network solution. 
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2.3.2 Structure Preserving Energy Function (SPEF) 

For classical model of generators, the following energy 
function has been defined [54] for the post-fault network ; 

5 

+ 2 W.. (2.4) 
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and the components of the potential energy, W 2 are 
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In the above expressions, the quantities with the subscript ’o' 
indicate their initial estimates. is the imaginary part 



33 


of j , the element of network admittance matrix [Y] (excluding 
machine reactances) . and are the COI referenced rotor 
angles of generator i and the voltage angles at bus i respecti- 
vely. It has been shown in [54] that the time derivative of W 
on the post-fault trajectory is zero. 

Comments 

1) It has been shown by Padiyar et al [70] that the last two 
components of the potential energy can be expressed as half 
the changes in the sum of reactive power losses in all elements 
of the network (including machine reactances) . Thusj 

*a4+'*'25 ° it=l ^ ^ '’1=1 ■ j=l 

where Qq^ = reactive power generation (at the internal bus) of 
generator i 

= reactive power load at bus j • 

L = total number of elements in the network including 
machine reactances 


in Eq. (2.11) can be calculated using the expression : 


^i 


-Ei^di 


( 2 . 12 ) 


2) The governor and prime mover system usually has little 
effect on stability. However, if the prime mover system is 
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included, mechanical power will not be constant. This can be 
included in the structure preserving energy function by modi- 
fying V^ 2 jl to 

m t d©. 

- « ( 2 . 13 ) 

Just as with non-constant active power loads, the above term 
can be numerically integrated. 

It is shown in APPENDIX C that the time derivative of 
the energy function VV including the non-constant mechanical 
input term is zero. 


3) '^ 22 * component of potential energy due to active power 

loads, is in general path dependent. Integrating, by parts, 
^22 resolved into two components, i) path independent 

and ii) path dependent, as follows s 


'»22 = ,1, [fpi(Vi) - { <^1 fpi (vpat 

^O ^0 


J, [(fpi(Vi) ’io’ - { ''i 


(2.14) 


It is observed that the path dependent term is zero for 
constant active power load characteristics. For simplicity, 
the path dependent term can be neglected even for the general 
case. With this approximation, the energy function is modified 



35 


and can be used for prediction of region of stability. It is 
shown in [71 ] that the approximation is acceptable and give 
accurate results. 


2.3.3 Computation of SPEF 

The total energy and its kinetic and potential energy 
components are calculated using a single subroutine TEF. The 
various input variables to the subroutine are indicated in 
Fig. 2.3. and 6 , obtained as output of the subroutine 

INT , are synchronous referenced variables. The corresponding 
COI referenced variables are obtained using the relations : 

0. = <5. "“6 « ad. — “Od f 3.~1»2j ...^m (2.15) 

1 10*1 1 o 


Here 6^ and u)^ are respectively the position and the speed 
of COI defined as 
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where m = no. of machines 

= moment of inertia of the ith machine 

The voltage Ej. behind transient reactance is obtained as 
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Fia2.3 RELATIONSHIP OF SUBROUTINE VARIABLES 
FOR CALCULATING SPEF 



37 


It is to be noted that the classical model generator can be 
obtained by choosing large values of and Also is 

made equal to to eliminate transient saliency. 

The bus voltages and their angles to be used in the compu- 
tation of SPEF are obtained as output of subroutine PFSOL which 
solves the algebraic equations of the network. COI referenced 
bus angles are calculated from the synchronous referenced out- 
put angles using the relation : 

5q » 1~1»^» »»*»n (2. 18 ) 

where is the bus angle with reference to a synchronous 
reference frame (calculated in PFSOL). 

A flow-chart of the overall program for direct transient 
stability evaluation using SPEF is shown in Fig. 2.4. 

2.4 DIRECT EVALUATION OF TRANSIENT STABILITY 

Whenever a fault occurs in a power system, the total 
energy of the post-fault system, whose stability is being 
examined, increases with respect to its initial value (which is 
customarily made zero) . If the fault is cleared after some 
time, the total energy is non-increasing. The key idea of the 
direct method is that the transient stability can, for a given 
contingency, be determined directly by comparing the total 
system energy which is gained during the fault-on period with a 
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FIG.2.4 OVERALL PROGAM FLOW DIAGRAM 
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certain 'critical energy'. The potential energy evaluated 
at the unstable, equilibrium point (uep) closest (in terms of 
energy) to the post-fault stable equilibrium point (sep) [21] 
gives conservative estimate of the critical energy and the 
region of stability. Recently, the concept of controlling uep 
in determining the critical energy has removed much of the 
conservativeness in the results [30]. Computational simpli- 
fications are possible with the use of PEBS method [29,30] in 
evaluating the critical energy (see APPENDIX D for a descrip- 
tion of the PEBS method) . 

PEBS method is particularly attractive as it avoids the 
need to evaluate the uep's.. In this method, it is assumed 
that the critically cleared trajectory goes near the uep 
corresponding to the fault location. PEBS is defined as the 
surface formed in the angle space by the points corresponding 
to the first maxima of the potential energy (with respect to 
sep) where its time derivative is zero. Inside this region, 
the potential energy increases and in the immediate vicinity 
outside this region, it decreases. Thus, the instant at which 
the time derivative of the potential energy function changes 
its sign from positive to negative, can be interpreted to be 
the instant when the trajectory crosses the PEBS [29,30]. The 
value of the potential energy at this time is the critical 


energy. 
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2.5 NUMERICAL EXAMPLES 
2.5.1 Description 

Three realistic power system examples have been considered 
for testing the transient stability program. The first system 
considered is a 7-machine (CIGRE) test system [72] . The 
second system is the 10-machine (NEW ENGLAND) test system [73]. 
The third system investigated is a 13-machine (UPSEB) system 
[74], The single line diagram of the systems are given in 
Figs. E.l , F.l and G.l alongwith system and operating 
data in APPENDICES E, F and G respectively. 

It is assumed that a three phase fault occurs at the 
terminal of a generator, the fault is cleared and the line is 
instantaneously reclosed. Thus, the prefault and the post- 
fault system configurations are assumed to be the same. Two 
different load characteristics are considered. For simplicity, 
loads at all the buses are assumed to have identical characteri- 
stics. Critical energy is determined using PEBS technique. 

The following cases are considered : 

Example 1 : 7-machine 10-bus (CIGRE) test system 

a) Fault at bus i; b) Fault at bus 2; c) Fault at bus 3’, 
d) Fault at bus 5; e) Fault at bus 6; f) Fault at bus 7. 

Types of load characteristics : 

(i) + jk^V^ ; (ii) + 3^2^^ 
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Example 2 : 10~machine 39-bus (NEW ENGLAND) test system 

a) Fault at bus 1 ; b) Fault at bus 3 ; c) Fault at bus 4 ; 

d) Fault at bus 5 ; e) Fault at bus 6 ; f) Fault at bus 7 $ 

g) Fault at bus 8 ; h) Fault at bus 9 j i) Fault at bus 10 * 

Types of load characteristics : 

i) + jk^V^ ; ii) + ^k^V^ 

Example 3 : 13-machine 71-bus (UPSEB) system 

a) Fault at bus 1 ; b) Fault at bus 3 ; c) Fault at bus 6 i 

d) Fault at bus 8 j e) Fault at bus 15 ; f) Fault at bus 27 ; 

g) Fault at bus 29 ; h) Fault at bus 39 ; i) Fault at bus 44 ; 

j) Fault at bus 48 ; k) Fault at bus 64 ; 1) Fault at bus 66 ; 

m) Fault at bus 68 . 

Types of load characteristics : 

(i) kj^V^ + jk^V^ ; (ii) + jk^^ 

2.5.2 Results 

In all the cases, critical clearing time, t^^ is obtained 
both by prediction and digital simulation and the results are 
con^ared. This is shown in Tables 2.1 - 2.3 for 7-machine, 

10-machine and 13-machine systems respectively for load types 
i) and ii) alongwith the corresponding value of the critical 
energy, 



Table 2.1 Critical Energy, and Critical Clearing Time, tor 7-inachine system 


43 


c 

o 




P 


o 

iD 


o 




P P 


CM 

o 


o 


CM 

>» 

CD CD 





vO 


> 

JQ 

P P 


• 



ft 


CM 

T3 

P D 


O 

o * 


O ?k 

M 

0> B 


1 

1 

1 


•r-) 

CD 

P P 


UD 

o 

1 1 

lO 

1 


C 

Q iO 


P 

o 


CON 

4* 

•H 






lO 



CO 


1^ 

• 

« 


ft 


rH 

p 



o 

o 


o 



XI 









o 







j 

• o 








1 

r— 1 

o 


TO 






CD 

CD 

C 

0 







CO 

o 

CO 






O 


>-P TO 

O 






S 

M 

COP O 


o 

Tl- 


CO 



CJ 

u o s: 

o m 

o 

C> >(< 

5!^ 

CON j 

T5 

P 

CD C P 

M m 


00 


iT) 


CD 


c: :3 CD 

act. 

k 

» 

1 1 

ft 

1 

O 


w u. s 


o 

o 


o 

h4 

P 









<0 



p 

lO 


lO 



o : 



CM 

CM 


CO 



•H 


O 

:li(: 


00 >k 


p ^ — 

•H CD M Z 


P 

« 

o 

ft 

t 

CO 

ft 

1 


U { 

: o ^ 


CM 

p 

CO 


o m 


P 

p 


rM 




c 








> 

o 


o 

lO o o 

in 

iD 


JD 

* ; p 


o 

lO O t1* 

lO 

uf) 


T5 

^ p p 

CO CO 


• 

k • 9 

nO 

ft 

00 

ft 

CM 

CD 

r-' P p 


o 

o o o 

Q 

O 


C 

*. p :o 


1 

1 1 r 

1 

I 

•o 

•H 

‘ CP s 


lO 

O lO iD 

O 

o 


CO 

>■-. *H ‘H 


o 

lO 00 00 

lO 

in 

+ 

P 

Q CO 


CO 

^ Sl“ ^ 

vO 

00 


P 




k 

'9 ft 

« 

ft 

CM 

O 



o 

o o o 

o 

o 

fH 










u 


TO 







CD 

C 

0 






«« 

CO 

o 

<0 








>••■-} "D 

o 






rH 

u 

CJ>-P O 


o 

vo lO r- 

ON 

o 

CD 

u 

M O X 

o u. 

C30 

^ 00 00 



TJ 

p 

<U C 

UJ 

00 

^ ^ 

vO 

00 

O 


C 3 <U 

acu 

k 

• 

ft '• 

ft 

ft 

s 




o 

o o o 

o 

o 

x> 

p 








CO 

CO 








o 

O 

>» 


r- 

§8 9 ! 

CM 

o 

*-4 

IP D> 


lO 

P 

lO 



-1 


CM 

ov vO CO 

CM 

CP 


P CD H :D 


00 

00 o 

ON 

CO 


IH 

Z OCL 


• 

ft 

ft ft 

ft 

ft 


O rjJ ^ sZ; 


lO 

IT) OD o 


in 





p 

p p p 

P 


C 









CD < 

3 








j CD 



p 

<N O') lO 

vD 

i> 

« 

1 









1 

3 








1 •« 

4 








p p 



p 

2 

3 

5 

vO 

0- 

P CO 








D O 


i 

CO 

CO cn CO 

CO 

CO 

1 CO 

0 



o 

D : 

3 :3 

D 

D 

1 U^.-4 



GO 

CQ CO CQ 

CQ 

CQ 


0) 

a> 


> 
c 
o 
, o 


4^ 

O 

c 



(f> 

CD 

o 

xi 



M 

O 

•P 

OJ 

)i 


* 



Table 2.2 Critical Energy, W and Critical Clearing time, for 10-machine System 
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The swing curves, the variation of the total energy and 
its components and that of the components of the potential 
energy function for the stable and unstable conditions are 
presented as follows : 

Figs. 2.5 to 2.10 For example 1-d) (i) 

Figs. 2.11 to 2.28 For examples 2-c) (i), 2-g) (i) and 

2-h) (i) 

2.5.3 Discussion 

It is observed that the predicted critical clearing time 
from structure preserving energy function is in good agreement 
with that obtained from digital simulation. It is further 
observed from Table 2.3 that using the modified energy fun- 
ction where the path dependent terms are neglected, gives more 
accurate prediction of critical clearing time. It is to be 
noted that the path dependent terms are identically equal to 
zero for constant active power loads. The results given in 
Table 2.3 are based only on the modified energy function. 

It is observed from the Tables that in all the cases 
considered for the three systems, the critical clearing time 
(both by prediction and digital simulation) is less for the 
constant active power loads compared to constant impedance 
type loads. This observation tallies with that made in 
[71] and indicates the importance of considering the voltage 
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dependent loads. In the case of nonlinear loads, there 
could be problems of convergence in network solution for 
some cases, indicating that the solution may not exist. 

These cases are indicated in the Tables. The problems of 
convergence with nonlinear loads can be avoided by dynamic 
representation of loads which is more close to the reality. 
However, this aspect is not investigated in this thesis. 

From the swing curves it is possible to predict the 
mode of instability. Invariably, the generator at the ter- 
minals of which the fault occurs, accelerates and goes out- 
of-step. In addition, there are other generators which may 
also accelerate and some generators which decelerate. In 
general, it is not possible to predict, in advance, which of 
the remaining generators (excluding the faulted one) .accele- 
rate or decelerate. 

In all the cases considered, the total energy increases 
during the fault-on period and remains constant after the 
fault is cleared. This is as expected and confirms the vali- 
dity of the results. It is interesting to observe that the 
system stability can be predicted from the variation of poten- 
tial and kinetic energy conponents. After the clearing of the 
fault, the kinetic energy and the potential energy oscillate 
and the magnitude of these oscillations are limited when the 
system is stable. When the system is unstable, the kinetic 
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FIG. 2.17 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 



FIG.2.18 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 
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Foult at bus 8 



FIG.2.21 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 
FOR 10-MACHlNE SYSTEM . UNSTABLE CASE 



FIG.2.22 VARIATION OF TOTAL ENERGY AND ITS COMPONENT: 
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Fault at bus ; 9 



FIG.2.25 VARIATION OF TOTAL ENERGY AND ITS COMPONEN 
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energy has a component that increases with time, while the 
potential energy has a component that decreases with time. 
Examination of the various components of potential energy 
shows that V}2^ or W22 or both can diverge. is the energy 

supplied from the generator rotors, while ^22 energy 

consumed in the loads. The difference between the two ener- 
gies is responsible for the increase in kinetic energy. It 
is also observed that ^^21 is negative if the generators tend 
to accelerate and positive if they tend to decelerate. 


W23, the energy cort¥)onent due to the reactive power in 
the loads, is not appreciable compared to other components. 
However, '^2^ + W25 which corresponds to the energy stored in 
the machine reactances and transmission lines is a signifi- 
cant component of the potential energy. This component is 
always non-negative and oscillates in response to the swings 
of the generator rotors. It is interesting to observe that 
all the components of the potential energy go through the 
minimum at the same instant of time when the system is stable. 
At these instants of time, the rotor angles are close to their 
equilibrium values. It is also interesting to observe that the 
frequency of oscillations in the magnetic energy (W24 + ^25) 
increases when the system becomes unstable. Also, for a large 
system, the peak magnitude of this energy for the unstable 
case is usually higher than that for the stable case. 
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2.6 CONCLUSIONS 

In this chapter, an existing dynamic simulation program 
[66] is modified to include structure preserving model of the 
system enabling detailed load representation in transient 
stability analysis. The program is then augmented by the 
addition of two hew subroutines to incorporate the computation 
of SPEF. The resultant program can, thus, be conveniently 
employed to evaluate the transient stability of power systems 
by i) digital simulation method and ii) direct methods using 
the SPEF. The transient stability program is tested on three 
realistic power systems. The following conclusions are drawn : 

1) Program output during the fault-on period can be analysed 
to determine the critical energy, and the critical 
clearing time t^^* 

2) The predicted value of t__ using SPEF agrees closely with 

V X s, 

that obtained by digital simulation in, practically, all 
the cases considered. 

3) The variation of energy components can give significant 
information about the system behaviour. In particular, 
the stability can be inferred from this instead of swing 


curves. 



CHAPTER 3 


A STRUCTURE PRESERVING ENERGY FUNCTION WITH KNOWN 

MODES OF INSTABILITY 


3.1 INTRODUCTION 

In Chapter 2, the application of structure preserving 
energy function to predict the critical clearing time was 
presented. This was based upon the computation of critical 
energy using PEBS method. The advantage of using PEBS method 
was that it avoids the calculation of controlling uep which 
can be quite complex [32] . 

The application of PEBS is based upon the assumption 
that the system is stable if the kinetic energy at the time of 
clearing gets converted into potential energy. However, in 
multi -machine systems, it is observed that the total kinetic 
energy need not be fully converted into potential energy [59 ] 
This is because of the relative motion among the generators 
that remain in synchronism. An iterative PEBS method has 
been suggested to account for this effect by modifying the 
critical energy or subtracting the conponent of the kinetic 
energy from the total kinetic energy at clearing, that does not 
lead to system separation [73]. The kinetic energy correction 
is equal to the minimum of the kinetic energy after the fault 
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is critically cleared. Implementing this correction, adds 
to the computational burden in predicting the region of 
stability. 

It has been suggested [32,59]that the kinetic energy 
responsible for system separation can be obtained by considering 
the mode of instability. For example, a three phase fault at 
the generator terminals can lead to the separation of the 
faulted generator from the rest of the generators which remain 
in synchronism. In this case, it is the kinetic energy com- 
ponent corresponding to the relative velocity between the 
generator going out-of-step and the remaining generators, that 
is contributing to system separation. This is an intutive 
concept that has been utilised in providing the kinetic energy 
correction in [59]. In [32] a two machine equivalent has been 
proposed to analyse this mode of instability. 

An exact derivation of the energy function based on mode 
of instability has not been presented so far. In this chapter, 
a novel structure preserving energy function is derived assum- 
ing that instability occurs due to the separation of the gene- 
rators into r coherent groups where, r is greater than or equal 
to 2. However, in practice, it is not possible to know the 
mode of instability apriori, yet if the fault is assumed to be 
at the terminals of a generator, then it is reasonable to con- 
clude that the faulted generator separates from the rest at the 
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onset of instability* This is a valid assunption in predicting 
•first swing’ stability using PEBS method. 

The application of the new SPEF for direct stability eva- 
luation to three realistic system examples is also presented. 

3.2 SYSTEM MODEL 
3.2.1 Assumptions 

In transient stability evaluation by direct methods using 
SPEF, the following assumptions are usually made : 

i) The generator is represented by a constant voltage behind 
the direct axis transient reactance (constant flux 
linkage) . 

ii) The governor action is neglected which implies that the 
power output from generator and turbine remains essen- 
tially constant during the transient. 

iii) Damping is ignored. 

iv) Transient saliency of the generators is neglected. 

v) Transmission lines are assumed to be lossless. This is 
generally true as the extra high voltage lines in a power 
system have high X/r ratio. 

The above assumptions are useful in the study of transient 
stability for only the 'first swing’ or for periods of the 
order of one second [ 2 ] and the error introduced by them are 
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small. However, the above assumptions can be relaxed without 
complicating the analysis. Assumption (v) has been relaxed 
in Chapter 5 for the development of system model and energy 
function to consider the transmission line resistances in 
stability study, while assumptions (i) and (iv) are relaxed in 
Chapter 6 for detailed modelling of machines. 

3.2,2 Generator Model 


3, 2. 2.1 General 

Consider a ’m’ machine system supplying 'n* nonlinear 
voltage dependent loads as shown in Fig. 3.1. Under the assump' 
tions mentioned in Section 3.2.1, the motion of the ith machine 
can be described by the following differential equations : 



0) 


syn 




= Pm-5“P^-5 

mi ei 


(3.1) 


E.V. sin(6.-.^.) 

where. P . — — ^ 5 i = 1,2, ...» m (3.2) 

Y ’ 

^di 

Eq. (3.1) is written with respect to an arbitrary synchronous 
reference frame . The coordinates of the rotational elements 
relative to COI of the system (global COl) describes exactly 
the synchronous behaviour of the system [36 ,37] .Global COI is 
defined as the inertia weighted average of all the rotor angles, 
1 .e . , 






T ransmission 
Network 


m 


m* 1 


m»2 


m*n 


Non-linear 
voltage dependent 
toads 


FIG. 3.1 A SYSTEM OF m MACHINES SUPPLYING POWER 

TO n NON-LINEAR VOLTAGE DEPENDENT LOADS. 
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1 m 1 ro 

^ “t = 4 “l <3.3) 

Following the formulation procedure outlined in ref. [so] » we 
get, 

•S- = 0)^ 
o o 

m 

“t “o = = ‘■cOI <3.4) 

The dynamics of the global CX)I reference are governed by the 
equations (3.4). By defining new angles and speeds relative to 
this reference 




(3*5) 


the system equations of motions become 


9. = w. 
1 1 




^i 

^mi “ ^ei *" ^C0I» 


i. ^ X f2. f * • • f ni 


(3.6) 


In view of the definitions in (3.3) and (3.5), the COI variables 

satisfy the following relations : 

m m 

i M* 9.! ~ 2 ” 0 

1=1 ^ ^ 1=1 ^ ^ 


(3.7) 
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3 .2.2.2 Aggregation of Generators in Coherent Groups 


When a disturbance occurs in the system, the generators 
in a coherent group swing together. Generally, the total 
number of coherent groups will be much less than m. In this 
case, it is possible to represent a coherent group by an 
equivalent generator and the order of the system dynamics is 
reduced. The equations describing a group a are derived 
bel ow : 


The local OOl of a group a can be defined as 



m 

2“ M. 9 . , 
i=l ^ ^ 




m 

M.to. , 
i=l ^ ^ 


m 2. = i:“ M. , 

T i=l ^ 


a 


where = number of machines in group a. 
Dynamics of group a can be described by 




(3.8) 


t a a 

o o 


,a •a 

U) 

T o 




a 


M~ co“ = S® (P . -P . ) 
^ '' 'mi ei^ 


i=l 


Ml ^COI» 


a 


1 , 2 , 


m • 9 


(3.9) 


Define new rotor angles and speeds of the machines in group a 
as 





i 


1,2, . . , , m^ 


(3.10) 


oc — \ fO. 9 


r . 
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Variables 6^ and 


satisfy 


the relations : 


"'a 

E = E M.wf « 0 , 

i=l ^ ^ i:=l ^ ^ 


OC ^ X 


(3.11) 


Substituting for u)^ in Eq. (3.6) from Eq. (3.10), the equation 
of motion of the ith machine in group a in the new variables 
*3^ and w? can be written as 

— Ot •CL 

“i“i = -“i“o + V - Pei - m; '’ool 


M, 


m 


a 


[ Z (P .-P .) •- 

*• \ * mot * AT / 


M 


a 


M.. 


‘‘isX ^ ~ ^rai’’^ei " ^CXff 


(3.12) 


substituting for w® from Eq. (3.9) 

Simplifying the R.H.S. of Eq. (3.12), finally we have, 

“l^ = Pmi-Pei - 3 P“cai - i .... 


M“ 


oc 1,2, ..., r 


(3.13) 


where. 


COI 


“a 

(Pml-Pei> - 


CC — i, f 2, f m • f X 


(3.14) 


Since the machines in group a are assumed to swing together, 
their angular rates are zero, i.e.,w®»0(i=l,2, ..., * 

Thus, Eq. (3.13) can be written as 


^i a 

^ “ ^mi'^ei “ j^a ^COI » ^ 

ot * 1 ,2 , . . . , r 


» “i, 


a 


(3.15) 
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Group, a is described by the dynamical equations (3.9) and the 
algebraic equations (3.15). It is to be noted that there are 

algebraic equations which can be used to 
solve for (m^-l) independent angle variables 0 “. It can be 
shown that the variables corresponding to the local CX)!* s of 
r groups satisfy the relations : 

r m 

2 M? 9 “ = S M. 9. = 0 « ^ (3.16) 

a«l ^ ® i=l ^ ^ 

« 

r m 

2 M- u = Z M. 0). = 0 (3.17) 

a=sl i=l 


Case of Two Groups 


This is a special case with r = 2. It is an important 
case as most of the times, at the onset of instability, one 
group (which may consist of a single generator with fault at 
its terminals) separates from the rest of the generators which 
may be clubbed into a single group. The dynamics of the 
COI's of the two groups are described by 





(P„.-P .) 

ini QX 



§ 






. ^^mi^^ei^ ‘ ^COI 

1 1 


(3.18) 


(3.19) 
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If the machines in groups 1 and 2 are assumed to swing toge- 
ther, the equations of motion of the machines in each group 
can be written as 


° ^mi~^ei “ j^i ^COI* ^ '*** '^1 


\ 2 


mi ei 


m: 


(3.20) 

(3.21) 


If group 1 consists of only a single machine (say, k), then 
Eqs. (3.18) and (3.19) reduce to 


“k 

^k^^k “ ^mk"^ek“ ^COI 
®T-k ^ ‘^T-k 


(3.22) 


^T-k “l-k ^ ^^mi'^ei^ " ^COI 

?^k 


(3.23) 


Here group 2 is labelled as (T-k) • Based on the assunption 
that the machines of this group swing together, their equa- 
tions of motion can be written as 


M, 

0 " Pm. P "" )\A 

mi ei M 


Z , ^COI,J-k’ ^ ...» m 

T-k ^ Ic 


(3.24) 


ra 


where P 


= 2 (P . -P . ) 

' vnn ' 

i=l 


COI , T- k . _i '■ ml ■ ei 


(3.25) 
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3.2.3 Load Model 


The real and the reactive power demands are modelled as 
arbitrary functions of respective bus voltage. Thus, 






(3.26) 

(3.27) 


3.2.4 Power Flow Equations 


At a node i corresponding to the terminals of a generator, 
the complex power injection can be written as (neglecting 
resistance) 


S. = 




= + (^, [ S t (-j B )]), i 


= + S2 


1,2,'.. *,m 

(3.28) 


The first term on the R.H.S. of Eq. (3.28) corres- 

ponds to power flow through the reactance of machine i conne- 
cted to node i and the second terra corresponds to the power 
flow in the lines connecting node i to node j (excludes the 


internal nodes of generators). In the above expression, 

B. . = Ira [Y. .] where [Y] is the admittance matrix of the net- 
ij i*3 

work (excluding machine reactances). 


The term can be written as 
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-j<p. -je. 

3 (p. V. e ^ - E. e ^ 
* V. e M 3 


E^V. e 

- j -j -Li- 


^di 


^di 


X 


di 


E.V. sin(<P.-S. ) V? 

-i-i + j [i 


E^V^ cos(<P^-e^) 


^di 


3 


x^. 

di 


X' 

^di 


(3.29) 


The term $2 can also be written as 


= 


j<fi. n -j<}) 

V. e ^ [ 2 V. e ^ (-jB..)] 
1 1=1, J 


n 




= j [ E B..V.V. e’ 3 

jrsl ij J 


n 


n 


« I B. -V-V.sin <P 
j=i ^ ^ 


ii - ^ "ij 


(3.30) 


where <P 


ij 



If there are no machines connected to node i, the complex 
power injection is given by 


s. = 7 . 1 ! 

1 11 


V. C 2 V* (-j B. . 
1 AzzX 


)] 


n 


n 


cos 


= Z B..V.V. sin <P. ,-j Z B..V.V u.., . . , 
jsirX J j=l J 

i = ra+1, ...» n (3.31) 
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Thus, in view of Eqs. (3.28) to (3.31), the real and the 
reactive power injections at bus i of a lossless system can 
be written as 

sin(^j^- 0 ^) n 


= 


"di 


£ B**V«V* sin r 1,2, •••, m 

(3.32a) 


n 


== £ B. -V. V. sin . , i = m+l , nH-2 , . . • , n 

ij 1 j ij 


V?' 


Qi = 


E.V. cos(©.-(p. ) n 

_i - Z cos<P.., i 

’'il 3=1 " 3 

n 

-SB. -V. V- cos <P . . j i = m+l , m+2, . ♦ . * , n 

4— 2^ xj 1 j ij 


(3.32b) 


1 ^ 29 .*«^in 
( 3.33a) 

(3.33b) 


Since the net real and reactive power at bus i is zero, we get 
from Eqs. (3.26) to (3.33) the following power flow equations : 


^i + ^pi^^i^ “ ^ ..., n 


^i ^qi^^i^ “ ^ ~ •••» n 


(3.34) 

(3.35) 


3.3 STRUCTURE PRESERVING ENERGY FUNCTION (SPEF) 

Consider the energy function,, defined for the post-fault 
system, as 


W = W2 = Wji 


(3.36) 


1 *** ot cc ^ 

where ^ii= kinetic, energy =5 ^ 


(3.37) 
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The conponents of the potential energy, W 2 are 
m 


W2i = 

- Z 
i=l 

(Si- 


n 

t 

^22 “ 

Z 

i-1 

•/' *’Li 



0 


n 

^i f 

*'^23 “ 

S 

i-1 

/ — 
^io 




dt 


gi 1 

X . 

1 


dx. 


W 24 = 


2^ [E?+V?-2E^Vi cos(e^-<P^) - (E^+V^^ 

cos(0i^“c^o^)3 


(3.38) 


(3.39) 


(3.40) 


(3.41) 


n n 


«25 = - 2 ‘'ij-VioVio ^3o> (3-42) 


The subscript ' o' in the above expressions indicates the quan- 
tities at the initial equilibrium. 

The time derivative of the energy function W in Eq. (3.36) 
can be shown to be zero as follows : 

Taking partial derivative of we get, 


— , a = l,2, ...,r (3.43) 

da)“ 

o 

Partial derivatives of W 2 w.r.t. ©^, and t are 


**' 2 - p 
6©^ ^mi 


E^Vj^ sin(e.-<PjL) 


X 


= -P 


di 


mi 


P • f 

ei^ 


i “ 1 >2 , 


* • f 


m 


(3.44) 
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dW2 


1 cos(0.-<P.) n 

T: ■ -L Bl/iVjCOs 

3 j- 


i “di 


di 


X X ^2 f « • • ^ m 


in 

~ V7 ■" ^ ~ nH-l,nH-2, n 


Vi ^^Li ■'' 


X. — 1 ^2 y • * • 9 n 


» 0, in view of Eqs. (3.27) and (3.35) 


(3.45) 


dW 2 sin(ei-<j).) n 

75 + 2 B. .V.V. sin V.., i = 1^2, m 

i di j=l ^ J 


n 


® 2 BijViVj sin ^i^, i = m+l, nH-2, ..«, n 

j “*2 


3= Pi* in view of Eqs. (3.32) 


(3.46) 


dW^ n 


d <P, 


ST 


= S P 


i^l "Li dt 


(3.47) 


From Eqs. (3.43) and (3.44), we get, 


dWii m dW 2 d 0 i r dWii da)“ r ®a ^v ^2 

^ i=l ^ ° a=l ^ «=1 1=1 ^ ^ 

r 

= S [m“ 0)“ w® + 2 (-P_ + P^.) ui.] 

a=l Too i^i mi ei' i-* 


1 
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m 


[M^ 0 )^ 0 )^ + + w“)], in view of 

Eq. (3.10) 


I' ^ "'a '^a 

2 ( [m“ w“ - Z (P .-P .)] a>“ + Z (-P .+P .) 

a=l T ° i=l ■'’^1 ei 0 ^ mx ex^ x' 


r M^“ 

■ a=l lil m“ '“I 


m ., — a 

a M.w. 

- E — :::: — ^rOT^ » view of Eqs. (3.9) 


and (3.15) 


= 0, on substituting from Eqs. (3.1l) and (3.17) 


(3.48) 


From Eqs. (3.46) and (3,47), we have, 


n dWo d^. 

< ^ — ar 

i=l d<P^ 


dWrt n d <P. 

’ ^ ^ ° i=l ^ 


= 0, in view of Eqs. (3.26) and (3.34) 
Hence, |jf = |e ("ll *2^ 


(3.49) 


dW^l^ m dW2 dS.. 


n dW^ dV^ 

“ IV 


= •^dT' + ^ 

1=1 1 1—1 1 

n dw„ dgj. dWo 

i"* JL ^ 


0, in view of Eqs. (3.45), (3.48) and (3.49) 

(3.50) 
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Comments 

1) The terms in Eqs. (3.38) to (3.42) of the potential 
energy function, W2 in Eq. (3.36) may be physically inter- 
preted in the following way : 

W22_ = change in potential energy due to mechanical input 
W22 ~ change in potential energy due to voltage dependent 
active power loads 

W23 = change in potential energy due to voltage dependent 
reactive power loads 

W24+W25 = change in magnetic energy stored in machine 
reactances and transmission lines. 


2] The kinetic energy term, Wjj_in Eq. (3.37) can be derived 
from that of the energy function developed by Padiyar et al. 
[54] as shown below : 
in 


r 

W. = w E M.w" = ^ 2 2 

^ ^ i=l ^ ^ ^ a=l i=l ^ 


1 .. .2 1 


(3.51) 


Substituting for in Eq. (3.51) from Eq. (3.10), we have. 


m 


ax2 


W, = i 2 E M. (w? + 0)“) 

^ ^ a«l i=l ^ ^ 

n\ 

= I f 1 “ + 5 S 2 


a=l i=l 
r ""a 

+ E 2 M.w? o)“ 

1 # «i jL mL Ki/ 

1=1 


a=l i=l 
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_ 1 


1 r , r r "“ct 

5 2 E + i E m“ (a)“)^ + E a)®( S M.w“) 

^ a=l i=l 1 ^ 2 T 0^ 1 1/ 


r a 


- S' E E M^Cok)^ + s E M~ (o)®)^, in view of Eq. (3.11) 

^ a=l i*l ^ ^ ^ a=l ^ ® 

(3.52) 


The first term on the R.H.S. of Eq. (3.52) corresponds 
to the total energy due to inter-generator swings in each of 
the r areas and does not contribute to the system separation. 
This portion of the kinetic energy need not be absorbed by 
the network for stability to be maintained and thus can be 
ignored. This results in an accurate estimation of the 
transient energy component responsible for system instability. 
Hence, 




(3.53) 


3) The kinetic energy term for the case of two groups can be 
written from Eq. (3.37) as 




(3.54) 


Using the relation in Eq. (3.7), W^^j^in Eq. (3.54) can be 


rewritten after some algebraic manipulation as 
11“ 2 "‘T 


Wn - ^ Mj (Uq) + j M-j. ( — 

M~ 

1 ,,1 ,1x2 /, ^ "‘Tx 

= ^ Mt (w ) (1 + ■^) 

Mj 

, Mi M? M_ wj o 
1 X,jr (.^)2 


5 “m: 


m: 


(3.55) 
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Similarly, it can also be shown that 


V/, 


1 2 
, m:: 

ii"" 5 ‘ ' 




U) 


0x2 


m; 


(3.56) 


The relative velocity of groups 1 and 2 is given by 
1 2 1 

~ view of Eq. (3.7) 

Wlj 




Also, 0 )^ 
o 



M 


M 


T ,2 

r% 


(3.57) 

(3.58) 


Utilising the relations in Eqs. (3.57) and (3.58), Wj^^^is 
obtained from Eq. (3.55) or (3.5.6) as 


1 Mi 1 o o 

*11= 5 t“o - 


(3.59) 


Similarly, when a single generator k separates from the rest 
of the generators belonging to group T-k the kinetic energy 
term can be given as 

1- q 

^ll~ ^‘^k ” ‘^T-k^ (3.60) 

4) For constant real power load, the eomponent W 22 

Eq. (3.39) of the potential energy function, W 2 can be directly 

computed from 
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^22 = 


n 

Z 

i«l 


- 


'»io) 


(3.61) 


In this case, the energy function is a Lyapunov function with 
the property that its Hessian is positive definite in the 
region of stability surrounding the stable equilibrium point. 

As discussed in Chapter 2, the expression given in Eq. 
(3.61) can be used as an approximation even for general case 
of nonlinear loads. However, the tinie derivative of the 
modified energy function is not equal to zero and is given by 


^mod 


n 

* Z 


i=l 


dV. 

''i W 


(3.62) 


It is observed that the modified energy function neglecting 
path dependent term gives quite accurate results [71]. 


5) It has been shown by Padiyar et al. [70] that the compo- 
nents W 24 and can be expressed as half the sum of rea- 
ctive power loss in each element of the network, i.e., 


, L , m n 

%4 ^25 * 2 ^3i “ 


(3.63) 


where 

* reactive power generation of machine i 



L = number of elements in the network 
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The R.H.S. of Eq. (3.63) is easily calculated at the end of 
power flow solution at each step during the transient, thereby 
simplifying the computation of the overall energy function. 

3.4 COMPUTATION OF STABILITY REGION 
3.4.1 General 

An index of the region of stability surrounding the 
stable equilibrium point ( sep) is the critical 
clearing time for a particular disturbance, typically a three 
phase fault at the terminals of a generator. The critical 
clearing time is predicted by confuting the critical energy. 

If at the time of clearing, the total energy is less than the 
critical energy, the system is stable and unstable otherwise. 

The critical energy for a given disturbance has been 
defined to be the energy evaluated at the controlling uep that 
lies on the stability boundary nearest to the sustained fault 
trajectory [30]. The evaluation of the controlling, uep can be 
avoided by the use of PEBS method (see APPENDIX D) . This 
approach is analogous to the equal area criterion used in 
two machine systems. In a two machine system, the potential 
energy reaches a maximum along the critically cleared traje- 
ctory at which time the kinetic energy (based on the relative 
velocity) becomes zero for conservative systems. Thus, it can 
be argued that the kinetic energy gained during the fault-on 



86 


period is completely converted into potential energy as the 
total energy remains constant for the post-fault system. 

The application of PEBS method to systems with more than 
two generators has following problems t 

i) It is not necessary that all the kinetic energy (with 
reference to COI as used in earlier formulations) gained 
during the fault-on period be converted into potential 
energy. Fouad and Stanton [59] has suggested subtracting 
the (first) minimum value of the kinetic energy along the 
critically cleared trajectory from the value of the 
energy at the time of clearing. However, this is a cum- 
bersome procedure and is not justified analytically* 

ii) While the potential energy maximum is same for the faulted 
trajectory and critically cleared trajectory, this does 
not extend to a multi— machine system. Since it is easier 
to evaluate the maximum potential energy along the faulted 
trajectory, an approximation is introduced by computing 
the critical energy from this maximum. The relative 
flatness of PEBS near the critically cleared trajectory 
is the justification for this, approximation [29] 

An iterative method of predicting the critical energy has 
to be used if this approximation is not valid [73]. 
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The mode of instability plays a major role in determining 
the amount of kinetic energy that is responsible for system 
separation. This is clearly brought out by the energy function 
derived in the previous section. While the potential energy 
is invariant and identical to that given in Chapter 2, the 
kinetic energy depends on the nature of coherent groups. 
Generally, there will be many coherent groups in a large system 
which may separate from each other. However, it may be con- 
ceived that, initially, one group of generators separate from 
the rest. As the time passes, these two groups may split into 
many subgroups. It is well known that energy methods as well 
as the equal area criterion are useful primarily in predicting 
stability in the first swing. While it is difficult to define 
* first swing' stability in a multi— machine system, it is rea- 
sonable to expect that the energy methods can predict the 
onset of instability where one or more generators start 
separating from the rest. 

It is to be noted that the dynamic trajectories of a 
multi-machine system are extremely complex and characterised 
by chaotic behaviour [19]. Thus, instability much after 

the fault clearing, cannot be ruled out. However, the presence 
of suitably designed controllers in the system to damp out the 
oscillations can overcome this problem. 
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3.4.2 Computational Procedure 

The evaluation of critical clearing time requires the 
simulation of the faulted system. Although the energy fun- 
ction is derived based on identifying the coherent groups 
and neglecting the relative motions within a group, simulation 
is performed without this assumption. Thus, the computer 
program, described in Chapter 2, is used with a minor modifi- 
cation in the computation of the kinetic energy component of 
SPEF. 

The most severe faults are the three phase faults at the 
terminals of generators. It is assumed that in the first 
instance, the generator at whose terminals the fault occurs, 
tends to separate from the rest. This is a reasonable assunp- 
tion and its validity is demonstrated by numerical examples. 

With the mode of instability thus assumed apriori, the 
corresponding kinetic energy component is used in the evalua- 
tion of critical energy. PEBS method is used with the usual 
approximation of conputing the critical energy from the faulted 
trajectory. 
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3 . 5 NUMERICAL EXAMPLES 
3.5.1 Description 

Three realistic power system examples considered in 
Chapter 2 are taken up again to illustrate the application of 
the proposed energy function to direct stability evaluation. 
Single-line diagrams of the systems are given in Figs.E.l,F.l and 
G.l respectively. The machine constants » line data and the 
operating data for the three systems are given in Tables 
E.1-E.3, F.1-F.5 and G. 1-0.4 respectively. 

It is assumed that a three phase fault occurs at the 
terminals of a generator. The fault is cleared followed by 
instantaneous reclosure of the line. The prefault and the 
post-fault network configurations are^ therefore, assumed to 
be the same. Two types of voltage dependent load characteri- 
stics are considered for investigation. For simplicity, load 
representation at all the buses are assumed to be identical. 

The following cases are studied : 

Example 1 : 7-machine 10-bus (CIGRE) test system 

a) Fault at bus 1 ; b) Fault at bus 2 ; c) Fault at bus 3 ; 
d) Fault at bus 5 ; e) Fault at bus 6 j f) Fault at bus 7.. 

Types of load characteristics : 

(i) + 3k2V^ 5 ^1 + ^*^2^ 



90 


Example 2 : 10-machine 39-bus (NEW ENGLAND) test system 

a) Fault at bus 1 j b) Fault at bus 3 } c) Fault at bus 4 ; 

d) Fault at bus 5 ; e) Fault at bus 6 } f) Fault at bus 7 ; 

g) Fault at bus 8 ; h) Fault at bus 9 ; i) Fault at bus 10. 

Types of load characteristics : 

(i) + jk^V^ ; (ii) k^ + Jk2V^ 

Example 3 j 13-machine 71-bus (UPSEB) system 

Fault at bus 3 (Gen .No. 2) ; 
Fault at bus 8 (Gen .No .4) ; 
Fault at bus 27(Gen.No.6) ; 
Fault at bus 39(Gen.No.8) ; 
Fault at bus 48(Gen.No.lO) i 
Fault at bus 66(Gen.No.l2) ; 


ki + jk2V^ 

3.5.2 Results 

Critical clearing times are predicted using the proposed 
SPEF and the SPEF given in Chapter 2 [54]. For all the cases, 
critical clearing times are also obtained by digital simulation 
and the results are compared. These are shown in Tables 3. 1-3 .4 


a) Fault at bus 1 (Gen. No.l) ; b) 

c) Fault at bus 6 (Gen. No. 3) ; d) 

e) Fault at bus 15(Gen. No. 5) ; f) 

g) Fault at bus 29(Gen. No. 7) ; h) 

i) Fault at bus 44(Gsn. No. 9) ; j) 

k) Fault at bus 64(Gen. No. 11); 1) 

m) Fault at bus 68(Gen. No. 13) . 

Types of load characteristics : 

(i) k^V^ + jk2V^ ; (ii) 
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Table 3.2 Effect of known Modes of Instability on the Predicted Critical Clearing 
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Table 3.3 : Effect of known Modes of Instability on the 
Predicted Critical Clearing Time, t^^ for 
10-machine System. 

Load Model ; i) 

Fault Gen. Critical t^^^Csec) obtained by 
Location No. Energy Energy Function Method Digital 

k (Modified SPEF) Simulation 

(pu) Neglecting Considering 

Mode of Mode of 
Instability Instability 


Bus 

1 

1 

9.2701 

0.262 

0.263 

0.260-0.265 

Bus 

3 

3 

9.2122 

0.232 

0.233 

0.230-0.235 

Bus 

4 

4 

9.8766 

0.214 

0.215 

0.210-0.215 

Bus 

5 

5 

3.3192 

0.169 

0.170 

0.175-0.180 

Bus 

6 

6 

10.5532 

0.238 

0.241 

0.240-0.245 

Bus 

7 

7 

9.5254 

0.225 

0.228 

0.235-0.240 

Bus 

8 

8 

10.1766 

0.228 

0.230 

0.230-0.235 

Bus 

9 

9 

2.3791 

0.108 

0.110 

0.120-0.125 

Bus 

10 

10 

23.5784 

0.628 

0.658 

0.570-0.575 
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for 7-machine, 10-machine and 13-machine systems respectively 
for load characteristics of types (i) and (ii). Tables 3.2 and 
3.3 present the results for 10-machine system for the following 
cases respectively ; 

i) When the path dependent term in the energy function is 
numerically integrated by trapezoidal rule 

ii) When the energy function is modified to approximate the 
path dependent term by a path independent term, as discussed 
earlier . 

As the results with the modified energy function are more 
accurate, only this case is considered for 13-machine system 
the results of which are reported in Table 3.4. 

Variables Q , V and (p are monitored corresponding to the 
critical energy for the contingencies investigated and load 
characteristic of type (i). These are presented in Tables 3.5 
and 3.6 for 7-machine and 10-machine systems respectively. 

Figs. 3. 2-3. 7 show the swing curves, the variation of 
total energy and its kinetic and potential energy components 
and also the components of the potential energy function, W 2 
for example 1-a) (i) for stable and unstable conditions. 

Swing curves for examples 2-a) (i) , 2-b) (i) ,2— d) (i) ,2-e) (i) » 
2-f)(i), 2-i)(i) and 2-a)(ii) for stable and unstable cases 

are presented in Figs. 3.8—3.21. The variation of total energy 
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Table 3.5 Variables Corresponding to the Critical Energy^ 

for 7-machine System (by Trajectory Simulation) 
Load Model : i) 



Q and ^ 

in radian 

and V in 

pu . 



Variables 


Fault at 

bus 




1 

2 

3 

5 

6 

7 


2.6147 

-0.2653 

-0.3709 

-0.1969 

-0.3966 

-0*3403 

Vl 

0.2716 

0.8231 

0.8081 

0.9254 

0.8568 

0.9745 

h 

-0.5229 

-0.4056 

-0 . 4866 

-0.3181 

-0.4896 

-0.4514 

®2 

-0.2774 

2.5939 

-0.1270 

-0.1616 

-0.4324 

-0.3549 

''2 

0.8319 

0.2705 

0.6234 

0.9211 

0.9098 

0.9905 


-0.4096 
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-0.4952 

-0.4416 

®3 

-0.3436 
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2.4998 
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-0.4508 
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''3 
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0.2301 
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0.8*994 

0*9655 

<^3 
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-0.3003 
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-0.4628 

®4 
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^4 

0.6413 
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0.9366 
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-0.4615 

-0.5340 

-0.3593 

-0.5248 

-0.4867 

®5 

-0.4972 
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-0.3637 

2.6787 

-0.6088 

-0.5172 

V 5 

0.8517 
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0.6702 

0.1611 

0.9130 

0.9800 


-0.5762 

-0.5027 

-0.5051 

-0.0926 

-0.6500 

-0.5832 
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1 

2 

3 

5 

6 

7 

®6 

-0.4910 

-0.5622 

-0.7059 

-0.3964 

2.4018 

-0.1860 


0.8107 

0.8558 

0.8195 

0.9331 

0.2351 

0.7844 

h 

-0.5700 

-0.5789 
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-0.4305 

0.0755 

-0.3579 


-0.7373 

-0.7841 

-0.9888 

-0.5079 

-0.0542 

2.4651 

^7 

0.9323 

0.9363 

0.8927 

0.9771 

0.7497 

0.3143 

0? 
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-0.7703 

-0.9545 
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-0.1490 

2.1653 

^8 

0.8432 

0.8519 

0.7711 

0.9278 

0.4966 

0.3960 

CD 

-0.6779 

-0.6793 

-0.8324 

-0.5009 

-0,2253 

-0.2686 

s 

0.7711 

0.6876 

0.3927 

0.8277 

0.7682 

0.8190 

05 

-0.5841 

-0.5190 

-0.6329 

-0.4147 

-0.5678 

-0.5313 

''lO 

0.7434 

0.7852 

0.6987 

0.5094 

0.8451 

0.9580 


-0.5643 

-0.5114 

-0.5480 

-0.3454 

-A. 6208 

-0.5646 



Table 3*6 Variables Correeponding to the Critical Energy, for 10-machine 
System (by Trajectory Simulation) 

Load Model : i) k.V^ + jk^V^ 
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its kinetic and potential energy components and also of the 
components of the potential energy function, for all these 
examples except examples 2-d)(i)and 2-i)(i) are shown in 
Figs« 3.22-3.41 for stable and unstable conditions. 


The variation of total energy and its kinetic energy 
component obtained using the proposed SPEF and the SPEF given 
in Chapter 2 are compared for stable and unstable cases and is 
presented as follows : 


Figs. 3.4 and 3.5 
Figs. 3.22 and 3.23 
Figs. 3.26 and 3.27 
Figs. 3.30 and 3.31 
Figs. 3*34 and 3.35 
Figs. 3.38 and 3.39 


for example l-a)(i) 
for example 2-a)(i) 
for example 2-b)(i) 
for example 2-e)Ci) 
for example 2-f)(i) 
for example 2-a)(ii) 


3.5*3 Discussion 

It is observed from Tables 3*1 - 3.4 that the considera- 
tion of the mode of instability in defining the kinetic energy 
responsible for system separation, results in improvement 
(although slight) in the predicted values of critical clearing 
time in, practically, all the cases. The notable exception is 
the 10-machine system for fault at bus 10 (example 2-i)(i)). 
However, it is interesting to note that machine 10 does not 


accelerate 
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FIG.3.3 SWING CURVES FOR 7-MACHlNE SYSTEM, UNSTABLE 
CASE 


Energy in pu 



V,I Vjo; 
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Fault at bus 1 

[*— = Fault cleared at 0.3 95 sec Load Model i) j K 2 V^ 




FIG.3.7 VARIATION OF THE COMPONENTS OF POTENTIAL ENERGY 
FOR 7-MACHINE SYSTEM, UNSTABLE CASE 



lie 




FIG.3.8 SWING CURVES FOR 10-MACHlNE SYSTEM^STABLE 
CASE 


Rotor angles (rad) — > , Rotor angles (rod) 
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FIG.3.10 SWING CURVES FOR 10-MACHlNE SYSTEM, ST ABLE CASE 
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FIG. 3.14 SWINS CURVES FOR 10-MACHINE SYSTEM, ST ABLE 
CASE 


Rotor ongles (rad ) — > , Rotor angles (rad) 
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FIG. 3.15 SWING CURVES FOR 10-MACHlNE SYSTEM, UNSTABLE 
CASE 














^'G.3.17 swing 

CASE 
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FOR 
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S'^STEM. UNSTABLE 
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FIG. 3.18 SWING CURVES FOR 10-MACHINE SYSTEM, ST ABLE CAS 


Fau(t Cleared 

at 0.575 sec Model : i ) k , v2 K2 v2 

Example ; 2-i)(i) 



®3 /~0 
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FIG. 3.19 SWING CURVES FOR 10-MACHINE SYSTEM^UNSTABLE 
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Rotor angles (rad) — > , Rotor angles (rad) 


Foult at bus 1 

Load Model -.ii) K'^* j K 2 V^ 

Example : 2- a) (i») 




FIG.3.21 SWING CURVES FOR 10-MACHlNE SYSTEM, 
UNSTABLE CASE 




FIG.3.23 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 
FOR 10- MACHINE SYSTEM .UNSTABLE CASE 





Energy in pu > Energy in pu 



FOR 10-MACHlNE SYSTEM, UNSTABLE CASE 


Energy in pu — > , , Energy in pu 




ENERGY FOR 10-MACHlNE SYSTEM, UNSTABLE CASE 



Energy in pu — > Energy in pu 



Energy in pu > Energy in pu 
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Fault at bus 6 

Load Model ; i) KiV^* j K 2 V ^ 



ENERGY FOR 10- MACHINE SYSTEM, STABLE CASE 



ENERGY FOR 10-MACHINE SYSTEM^ UNSTABLE CASE 



Fault at bus 7 



^ jr 


Fault at bus 7 



ENERGY FOR 10-MACHINE SYSTEM , STABLE CASE 



FIG.3.37 VARIATION OF THE COMPONENTS OF POTENTIAL 
ENERGY FOR 10-MACHINE SYSTEM, UNSTABLE 
CASE 



nrl 111 XRi^ii*3 
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lO.Or I pQult cleared at 0.230sec Fault at bus 1 

[ Load Model ; U) Ki ♦ j K 2 V^ 



_5.oLF1G.3.38 variation OF TOTAL ENERGY AND ITS COMPONENTS 
FOR IQ-MACHINE SYSTEM . STABLE CASE 


1 

10.0r ■— F oult cleared ot 0.235sec / 



[•I*j 


Fto. 3.39 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 
FOR 10-MACHINE SYSTEM.UNSTABLE CASE 



Energy in pu — > , Energy inpu 



FlG.a40 VARIATION OF THE COMPONENTS OF POTENTIAL ENERGY 
FOR 10-MACHlNE SYSTEM^ STABLE CASE 



FIG.3.41 VARIATION OF THE COMPONENTS OF POTENTIAL ENERGY 
FOR 10-MACHlNE SYSTEM, UN ST ABLE CASE 
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In general, the mode of instability can be quite complex 
as seen from the swing curves. Yet it would appear that this 
complexity does not affect much the predicted results. The 
major reason for this is that the correction in the kinetic 
energy term to account for the mode of instability is of 
secondary importance in comparison with the correction in the 
potential energy term. Accurate prediction of the mode of 
instability is essential when the critical energy is determined 
from the energy evaluated at the controlling uep as the con- 
trolling uep is related to the mode of instability. However, 
the accuracy of the PEBS method in determining the critical 
energy is not dependent on exact knowledge of the mode of 
instability. 

The total kinetic energy consists of two components; 
i) defined in Eq. (3.37), which is responsible for system 

separation and ii) ^^, 2 * component which corresponds to the 
relative motion of the rotors of the generators within any 
coherent group. While the energy function derived in . 

Section 3.3 assumed 't-o zero, this is not required as 

long as W 2^2 remains constant along the post-fault trajectory. 

As perhaps to be expected, ^^2 <^oes not remain constant (see 
Figs. 3.22, 3.26, 3.30, 3.34 and 3.38) and oscillates. This 
is a fact of crucial importance in explaining why the kinetic 
energy correction neglecting Wj ^2 adequate. 
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In stably cleared fault, it is expected that goes 
through zero following the clearing of the fault. Thus, 
the minimum of the kinetic energy, following the clearing 
of the fault, is equal to Wj^2* Fo^ad and Stanton [59] 
assume that the minimum kinetic energy is equal to W ^2 at the 
time of clearing. However, this is not correct as seen from 
Figs. 3.22, 3.26, 3.30, 3,34 and 3.38. It is observed that 
minimum kinetic energy is much larger than W ^2 at the time of 
clearing, even though passes through zero following the 
clearing of the fault in the stable case. This shows that 
part of the kinetic energy component is transformed into 
the component which corresponds to the relative motion of 
the generators within a coherent group. Thus, the correction 
of the kinetic energy neglecting Wj ^2 is not adequate. The 
kinetic energy correction by considering the minimum kinetic 
energy (following stably cleared fault) is more accurate. This 
is shown in Fig. 3.42 for the stable case for example 2-a)(i)» 
where the kinetic energy correction by the later method gives 
a critical clearing time of 0.263 sec, whereas the predicted 
time neglecting 0.249 sec (Table 3.2). 

It is interesting to note that for large systems (both 
10-machine and 13-machine) neglecting the path dependent term 
in the energy function (duo to the non-constancy of active power 
loads) giyes a very good result which makes the kinetic energy 
correction superfluous. Vifhile it is difficult to explain the 
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FIG. 3.42 KINETIC ENERGY CORRECTION FOR 10-MACHINE SYSTEM , STABLE 
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phenomena, it. simplifies the computation of the energy fun- 
ction. 

The examination of the energy function and its various 
components brings out the following points and supports the 
conclusions in Chapter 2. 

i) the instability can be judged from the variation of 

or W2« While Wj^ continues to increase, W2 continues to 
decrease (as the total energy is constant along the post- 
fault trajectory). Instability may occur after the first 
swing. 

ii) W23, the component of the potential energy corresponding 

to the reactive power loads, is insignificant compared to 
^22 active power loads. This shows that in- 

accuracies in modelling the reactive power component of 
loads have negligible effect in predicting the stability 
region. However, this cannot be said about the modellinc 
of the active power loads. 

iii) The potential energy terms and W25 are always non- 
negative and bounded even for unstable cases, although 
higher peaks are reached if the system is unstable. 

iv) The potential energy texm ^21 mechanical input, 

continues to decrease (or increase) monotonically when t\ 

system is unstable. 
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It is interesting to note that the oscillations in 
continue to decrease for the two cases (examples 2-e)(i) and 
2-f)(i)) even when the system is unstable. This is due to the 
fact that, in both the cases, all the generators tend to acce- 
lerate together and separate from generator 2 which decelerates. 
This implies that the mode of instability assumed (that the 
faulted generator separates from the rest) is not correct. How- 
ever, the effect of this error in predicting t - insigni- 

cr 

ficant. 

3.6 CONCLUSIONS 

In this chapter, a system model is developed based on the 
aggregation of generators into r coherent groups on the occu- 
rrence of a fault, each group having one or more generators in 
it. A new structure preserving energy function is then derived 

I 

assuming that the mode of instability is known apriori. Appli- 
cation of this energy function for predicting the stability 
regions for three sample systems shows that 

i) the kinetic energy correction neglecting compo- 

nent corresponding to the relative motion of the generator 
rotors within the coherent groups) gives slight improve- 
ment in the predicted value of critical clearing time. 
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ii) although the component of kinetic energy that leads 

I 

to system separation, goes through zero for a stably 
clear fault, it does not get converted completely into 
potential energy; part of gets converted into Wj^2* 

The exchange in these two components of energy, and 
Wi2> complicates the exact determination of stability 
region for large multi -machine systems. 

iii) neglecting the path dependent component of the energy 
function, due to the non-constancy of active power loads, 
appears to improve the accuracy of prediction even with- 
out kinetic energy correction. 

iv) PEBS method using SPEF appears to be quite robust in 
predicting the critical energy. The kinetic energy 
correction has, basically, a minor effect on the accuracy. 



CHAPTER 4 


DYNAMIC SECURITY ASSESSMENT USING STRUCTURE PRESERVING 

ENERGY FUNCTIONS 


4.1 INTRODUCTION 

While some attempts have been made to apply direct methods , 
based on energy functions on reduced network models, to on-line 
dynamic security assessment [64, 65], no attempt has been tried 
with structure preserving energy functions. The objective of 
this chapter is to present algorithms and results in applying 
SPEF for on-line dynamic security assessment. 

The security assessment is an important function in modern 
energy control centres to check whether the present operational 
state is normal. While the methods of static security assess- 
ment are well advanced and implemented in many control centres , 
the problem of dynamic security assessment is still not well 
defined. One of the objectives can be the checking of transient 
stability following a contingency such as a three phase fault 
at the terminals of a generator. 

The use of transient energy margin as an index helps not 
only in checking the stability, but also to rank the contingen- 
cies according to their severities. The computation of tran- 
sient energy margin requires the computation of critical energy 



141 


and the energy at clearing. The later usually requires the 
simulation of the faulted trajectory. Approximations of the 
faulted trajectory has been used to simplify the computations 
[64]. The critical energy can be determined using controlling 
uep or PEBS method. Both require again the simulation of the 
faulted trajectory. The use of SPEF, although attractive from 
the view point of preserving the structure of the system, 
requires more computational efforts in the simulation of the 
faulted trajectory. Both faulted and post-fault system net- 
works are solved at every step to compute the faulted tra- 
jectory and the energy function respectively. Efficient methods 
for speeding up the computation of transient energy margin 
would be required before SPEF can be seriously considered for 
on-line dynamic security assessment. 

The method proposed in this chapter assumes that a faulted 

generator accelerates with respect to the remaining generators 

and its rotor angle increases monotonically during the fault-on 

period. If it is assumed that the rest of the generators swing 

► 

together, then their rotor angles along with the bus voltage 
magnitudes and angles become implicit functions of the rotor 
angle of the faulted generator. The potential energy component 
of the structure preserving energy function also reduces to a 
function of a single independent variable. With this assumption 
the estimation of the faulted trajectory is avoided. 



The method is illustrated by considering the three system 

t 

examples of Chapters 2 and 3. 

4.2 TRANSIENT STABILITY MEASURES 

4.2.1 General 

The energy margin can be used as an index for assessing 
the robustness of a power system in a given configuration and/or 
predicting its transient behaviour under abnormal conditions. 

To the system operator, it offers a means of performing on-line 
dynamic security assessment. 

4.2.2 Transient Energy Margin (TEM) [61 » 63-65] 

The transient energy margin is defined as the difference 
between the critical energy and the energy at the clearing time. 
The index serves to indicate the severity of the disturbance. 

Thus, if is the value of the energy function, W at the 
clearing time and the critical energy, the transient energy 

C i 

margin, A vi/ can be defined [61, 63-65] as 

( 4 . 1 ) 

The transient energy margin, AVI indicates how close to insta- 
bility a system might be, for a given system condition and for a 
particular disturbance. If A W is positive, then the system is 
stable. If AW is negative, the system is unstable. 
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Normalised Energy Margin [32] 


An alternate measure of transient stability is the ratio 
of potential energy margin of the system to the kinetic energy 
at fault clearing. Denoting the kinetic and the potential 
energy components at clearing by and ^ 2 ^,^ respectively, 

the normalised energy margin, A may be defined [32] as 


AWn 


""cr ' ^^2cl 


W 


Icl 


(4.2) 


This is used to evaluate the system stability in the following 
manner : 

i) a value of AW^ greater than 1 indicates stability; if A 
is less than 1, the system is unstable. 

ii) a disturbance which gives the lowest value of consti- 

tutes the most severe disturbance. 


4.3 DYNAMC SECURITY ASSESSMENT 
4.3.1 General 

Dynamic security evaluation of power systems through the 
application of transient energy margin or the normalised energy 
margin, involves the computation of i) critical energy and ii) 

transient energy and its kinetic and potential energy components 

* 

at the instant of clearing the fault. Critical energy can be 
evaluated using the PEBS method which requires the simulation 



of the faulted trajectory. Although the application of this 
method avoids the computation of uep' s for evaluating the cri— 
tical energy » the computation of the faulted trajectory may 
constitute a limitation in its application to on-line dynamic 
security assessment. A method is proposed in this section 
which eliminates the need to determine the faulted trajectory 
and can speed up the computation of TEM and the normalised TEM. 

As discussed in Chapter 3, for a fault at its terminals, a 
generator (say, k) gets severely disturbed and initially tends 
to separate from the rest of the generators which may be assu- 
med to swing together. This implies that the generator k acce- 
lerates monotonically during the transient. Under the assump- 
tion that the remaining generators swing together, their motions 
can be described by algebraic equations. This permits the rotor 
angles of these generators and also the bus voltages and their 
angles to bo expressed implicitly as functions of the rotor 
angle of generator k. This also reduces the potential energy 
to a function of single independent variable corresponding to 
the rotor angle of the generator going out-of-step with respect 
to the rest of the generators. 

4*3.2 Mathematical Formulation 

Mathematically, the problem for computation of TEM can be 


formulated as follows 



i4i) 


Consider the modified energy function (neglecting the 
path dependent term) defined in Section 3.3, Chapter 3 [eq. 
(3.36)] as 

W = Wj^j^(ojj^) + W2(e., V, ej^) (4.3) 

where 0 is (m-1) vector of rotor angles (excluding the faulted 
generator) , V and 2. ai'e n vectors of bus voltage magnitudes and 
angles respectively. 

Under the assumption that all the generators except. the 

faulted generator swing together, it can be shown that the 

vectors 0 and £ are functions of Using PEBS method, the 

critical energy can be determined, if can be determined such 

* 

that W 2 ( 0 jj) is maximum along the faulted trajectory. 

The energy at the time of clearing can be determined if 
the values of 0^^ and at clearing are determined. Thus, 

W{t=tci) = (4.4) 

if 

4. 3. 2.1 Computation of and Critical Energy 

The dependence of vectors 0, V and £ on 0j^ is shown by 
the following equations which are taken from Chapter 3 : 

There is no loss of generality in assuming that k = m 
(this can be done by renumbering the buses, if necessary). 
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From Eqs. (3.24) which describe the motions of the genera- 
tors within the coherent group, we have, 

M. m 

9i = ° = Pmi - "el (®i- '»i) - f i)] 


(4.5) 


i ~ 1,2, ..., m— 2 


From Eqs. (3.34) and (3*35) which describe the power flow 
equations at n buses, we have, 

9jjj_2+i “ 0 “ ^i^^i’ ^i * ^i^ ^ ^pi^^i^ ^ ^ ~ 1»2, m (4.6a) 

^2m-2+i” Pi(V^,<P^) + fpi(V^), i = 1,2, . . . ,(n-m) (4.6b) 

^m-2+n+i” ^ ~ ^i^^i* ^i* ^i^ ^qi^^i^ ~ ly2, *.., m (4.7a) 

92m-2+n+l=° = ‘^1) + fqi(Vi) .1 = 1.2. ...(n-m) (4.7b) 


From Eq. (3.7), we have, 


g 


2n+m-l 


m 

= 0 = E M, 9 


i=l 


11 


(4.8) 


Eqs. (4.5) - (4.8) can be written in the compact form : 


g (©» y, V , e^) = 0 (4.9) 

In the above expression, g is a vector valued function of di- 
mension of (2n+m-l). The total number of unknowns 0., V and 2. 
are also 2n+m--l. Hence, given 0j^, we can solve for the unknown 
variables from Eqn. (4.9). 
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Denoting [S, V, SP]^ = and = u, the problem of 
determining the critical energy, W can be formulated as 

O Jl 

follows : 

Max W 2 (x,u) (4.10) 

subject to g(x, u) = 0 (4-11) 

Let the solution of the problem be u = u*. Then, 

^cr ^ (4.12) 

* 

v/here _x is the solution of 

g (x , u ) = 0 (4.13) 

Noto : In the presence of path depen.dent term in the energy 
function, numerical integration using trapezoidal rule 
is used to evaluate these terms and the energy function 
can be written as 

W(t=t^) = W^^(wJ) .+ (4.14) 

where superscript r indicates the value of the variable 
evaluated at t = t^ and similarly superscript r-1 _iiaidi- 
cates the value of the variable evaluated at t = • 

The method of computation of W is now modified as 

W Jm 

iollows : 

''''^cr “ -r-1’ ’^r^ 


(4.15) 
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subject to £ u^_^) 

= 0 

(4.16) 

g (Xpj 

= 0 

(4.17) 

and u_ 

r 

== Vl 

(4.18) 


where ~ ~ ’^r “ 

^ u^ can be arbitrarily chosen. 


4.3.2 .2 Computation of Energy at Clearing 

As mentioned earlier, this requires the knowledge of 
and a)j^ at clearing which can be obtained from 

M, m 

“^k^k ^mk ' ^ek ~ ^^mi”^ei^ 

M, 


m 

= - M;)("mk “ Pek) - m; , 1, ("’mi-Pei) 




= P 


ak 


(4.19) 


Under the assumption that the remaining generators swing toge- 
ther throughout the fault-on duration starting from t = 0 , 
Pak can be shown to be constant. Pj^^^ (i = 1»2, ..., m) is 
constant as classical model is assumed. Power output of 
generator k, = 0 during the fault. P^j^ can be calculated 
from 


Pel = Pel(®l> ^ “ 


(4.20) 
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where 


9 

i “ ®i " ®T- 

•k '^1 = f i - 8^.^ 


Define 

vectors 



-T 

B 

= [^1 Sj 

• • • 

(4.21) 


rH 

> 

II 

’•• ''m-1 ''m+l ••• 

(4.22) 

-T 
(P . 

= [»i 

••• Vl »iw-l ••• 

(4.23) 

It is 

to be noted 

that during the fault, = 0 and 

9^ = 0. 
m 


The vectors, V and <P can be obtained from the equations ; 

M. m 

h. = 0 = 2 

r**ic 1—1 

1 ” 1 , 2 , •••, ni "**2 

(4.24) 

^2+1 = ° = \'> * ^ “ 

(4.25a) 


*'2m-3+l = ° = ^ ° 

*’m-3+n+l= ° = Qi(\-Vi.^>+fqi(Vi), i = 1.2,...,m-l (4.26a) 
‘'2m-4+n+l=° = <3i‘Vi.\) + fqi(Vi) . 1 = 1,2, .. . ,(n-m) (4.26b) 


and finally, 



= 0 = 


(4.27) 


2n+m-3 


m 

E M, e . 
i=l ^ 1 


Eqs. (4.24) - (4.27) can be written in the compact form : 


ll (1» V.D = 0 


(4.28) 


where h has the dimension of (2n+m-3). Eqn. (4.28) can be 
solved to obtain ©, V and which are assumed to remain con- 
stant during the faulted period. Once these vectors are deter- 
mined, (i = 1,2, m-1) can be computed and this deter- 


“kcl ®kcl determined from 

m - ^ak ^cl 

‘^kcl " 


jccl 


®ko + 


’ak ^^cl^' 

2M,. 


where t^.^^ is the clearing time and ~ ®k = 0) . 


(4.29) 


(4.30) 


Substituting Eqs. (4.29) and (4.30) in Eq. (4.4), we can 
get the kinetic energy and the potential energy at the time of 
clearing. 


4.3.3 Computational Procedure 

The major computational effort is the determination of © 

to evaluate W . 0^, the rotor angle of the faulted generator, 

ex' K 


tv * 
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increases monotonically during the fault. The approach used 
here is to simulate the trajectory of 0 j^ during the fault 
and monitor the potential energy^ V ^2 until it crosses the peak. 
The trajectory of is simulated by starting from ®ko ^ 

©kCt = O) and increasing its value in small steps. 

The solution of the set of nonlinear equations (4.9) can 
be simplified by linearisation. 

Given g ^r-1^ ~ - (4.31) 


We can write, 


dg 


dg 


0 = £ (Xr» ^r-1^ ^ ^ ^"^r-1 ^ ^ 


'u=u 


r -1 


x=x 

u=u 


r -1 

r -1 


from which 

AUj, 


(4.32) 

(4.33) 


where [s] = [ [^] 

,d g d g 

The elements of and (^g^) are given in .^PPSdDIX H. 


Thus, given the solution of ^ at the previous step, the solu 
tion at the present step can be obtained as 


X = X . + 

— r — r-1 


AX, 


(4.35) 


where Ax^, is the solution of the linear Eq . (4.33). 



IJ'Z 


The accuracy of the solution of ^ is controlled by selecting 

the step size Au . 

r 

The algorithm for the determination of critical energy is 
as follows : 

Step 1 : Solve for x at t = 0 from 

2 =0 (4.36) 

using Newton's method. 

Step 2 : Compute W^Ct = 0) = ^o^ 

Step 3 : Set iteration count r = 1 and tiU = H 

r 

Step 4 ; Augment u^_^ hyAu^, i.e,, 

Uj.=Ur_i +iUj, 

Step 5 : Solve for A 2 <^ using Eq. (4.33), i.e,, 

= [S]AUj^ 

Step 6 : Update the value of 2£r_i» i.e«» 

X = X . + A x^ 

— r — r-l -r 

Step 7 ; Calculate 

AW^ = W2(t^) - W2( 

.If AW 2 is positive, set r = r+1 and go to Step 4. 
Otherwise, set stop. 
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4,4 NUMERICAL EXAMPLES 
4.4.1 Description 

The three system examples considered in Chapter 2 and 
Chapter 3 are again taken up to illustrate the application of 
the proposed method for dynamic security assessment of power 
systems. It is assumed that a three phase fault occurs at the 
terminals of generators which is cleared followed by instan- 
taneous reclosure of the line. This implies that the prefault 
and the post-fault networks are assumed to have identical 
configuration. For simplicity, all the load buses are con- 
sidered to have similar characteristics. The following examples 
are investigated : 

Example 1 s 7-machine 10-bus (CIGRE) test system 
Fault locations : 

a) bus 1 ; b) bus 2 ; c) bus 3 ; d) bus 5 * 
e) bus 6 ; f) bus 7 

Load models ; 

1) kj^V^ + jk^V^ ! ii) ‘‘l + 

Example 2 : 10-machine 39-bus (NEW ENGLAND) test system. 

Fault locations ; 

a) bus 1 5 b) bus 3 ; c) bus 4 ; d) bus 5 •, 
e) bus 6 f f ) bus 7 ; g) bus 8 ; h) bus 9 
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Load models : 

i) + jk^V^ ; ii) + jk^V^ 

Example 3 : 13-machine 71-bus (UPSEB) system 
Fault locations : 

a) bus 6 (Gen. no.3) ; b) bus 15 (Gen. no. 5)^ 

c) bus 27 (Gen. no,6) ; d) bus 39 (Gen. no. 8); 

e) bus 44 (Gen. no.9) ; f) bus 64 (Gen. no. 11) 

g) bus 68 (Gen. no. 13) 

Load models : 

i) kj^V^ + jk2V^ ; ii) k^^ + jk2V^ 

4.4.2 Results 


The critical value of the energy function is evaluated 
using PEBS method for all the three examples both by performing 
the simulation of the faulted trajectory and by the proposed 
method. The results are then compared and are presented in 
Tables 4. 1-4. 3 for examples 1,2 and 3 respectively. An 
arbitrary step size of H = 0.1 radian is used in the computa- 
tion of 0j^. 

The rotor angle and velocity of the faulted generator k.are 
obtained for example 2 using i) the proposed method and ii) 
trajectory estimation of the faulted network, for various fault 
locations and. load model (i). The results are compared and are 

given in Table 4.4. These results are obtained for an arbitra- 
rily selected clearing time, = 0.20 sec. 



Table 4.1 Comparison of critical energy, W for 7-machine system 

y w jL 
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Table 4.2 Comparison of Critical Energy, W._ for 10-machine System 
0 in radian and in p.u. 
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Table 4.3 Comparison of Critical Energy, for 13-machinG system 
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The kinetic energy and the potential energy components of 
transient energy at clearing are estimated corresponding to the 
clearing time, t^^^ = 0.20 sec using i) the proposed method and 
ii) the faulted trajectory simulation, considering load represen- 
tation of type (i) for example 2» The transient energy margin, 
AW and the normalised transient energy margin, A are then 
calculated and the contingencies are ranked based on their seve- 
rity. Results obtained using the proposed method are presented 
in Table 4.5 while Table 4.6 shows the results obtained by tra- 
jectory simulation. 

4.4*3 Discussion 

It can be observed from Tables 4. 1-4. 3 that the value of 
the critical energy obtained by the proposed method is generally 
in good agreement with that obtained by performing the trajectory 
simulation of the faulted system. For 7-machine system, the 
maximum and minimum differences in the value of critical energy 
are 3/» (fault at bus 3) and 0.39;^ (fault at bus 6) respectively 
(Table 4.1) with the loads modelled as (i)* same load 

model, the corresponding values for 10-machine system are 9 •58^ 
(fault at bus 4) and 0.59;^ (fault at bus 9) respectively 
(Table 4*2). For 13-machine system, the maximum and minimum 
values can be seen (Table 4.3) to be 3.11?^ (fault at bus 39) and 
0»92yi (fault at bus 64) respectively. The maximum and minimum 
differences in the value of the critical energy for load 



Table 4.5 AW and by the Proposed Method and Ranking of Three phase Faults 

for 10-machine System 

Clearing time, t , = 0.200 seconds 
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characteristics of type (ii) are less than those obtained fox 
the load representation of type (i) for the fault locations 
considered on all the three system examples investigated. The 
problem of convergence of the network solution is also encoun- 
tered with constant active power loads as in Chapter 3 and this 
is indicated in the tables. 

It can be seen from Tables 4.5 and 4.6 that the proposed 
method for computing TEM predicts accurately the stability of 
the system for all the disturbances considered. Although the 
values of TEM calculated by the proposed method do not exactly 
match with the values obtained from simulating the faulted tra- 
jectory, it is interesting to observe that the ranking of the 
contingencies given by both methods are identical. This shows 
the practical validity of the proposed method. 

Actually, the proposed method can be viewed as an exten- 
sion of the equal area criterion for multi-machine systems, 
based on the assumption that the faulted generator separates 
from the rest. Although a two-machine equivalent can be derived 
under this assumption for constant impedance type loads [32] , 
the procedure given in this chapter does not require such a 
derivation for simplifying the computation. Also, for stru- 
cture preserving models, it is not feasible to obtain two- 
machine equivalents, particularly for nonlinear type loads. 
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It is interesting to note that the accuracy of the 
proposed method is not critically dependent on the validity 
of the assumption about the mode of instability. For exanple, 
in the case of 10-machine system it is observed (see Chapter 3 
that the faulted generator accelerates along with other genera 
tors in the case of faults at buses 6 and 7. 

The major computational effort in the evaluation of TEM 

in the proposed method, lies in determination of W . Even 

cr 

here there is a significant saving in the computations than th 
usual method of integrating the faulted system equations. The 
use of SPEF requires the solution of the network equations 
twice at a step during the simulation of the faulted equations 
The proposed method requires the solution of the post-fault 
system equations and even this is simplified using linearisa- 
tion. In the numerical examples considered, it required abou 
20-25 steps to confute W in the proposed method, while it 
took 40-50 time steps to compute using the simulation of 
the faulted equations. Both faulted and post-fault networks w 
solved for each step during the simulation. 

Further computational savings are possible in the propose 
method if approximations can be introduced in the sensitivity 
matrices in each step. 

It is interesting to observe that the calculation of 
and in the proposed method do not require the knowledge o 
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faultsd ‘txa j ©c 'tory • Alsoj the coinputstion of energy at 
clearing requires the solution of one set of nonlinear alge- 
braic equations (to determine the accelerating power of the 
faulted generator) . 

4.5 CX)NCLUSIONS 

In this chapter, a novel method using SPEF developed based 
on an apriori knowledge of the mode of system separation, is 
presented to evaluate the transient energy margin. This method 
does not require the simulation of the faulted equations to 
compute critical energy and energy at clearing and is computa- 
tionally quite efficient. The application of the proposed 
method to three large scale sample systems yields encouraging 
results and lead to the following conclusions : 

i) the proposed method is the first step in the application of 
SPEF for on-line dynamic security assessment. The results 
obtained from this method is promising enough to continue 
the investigation for speeding up the computation required 
in the on-line dynamic security assessment. 

ii) the normalised transient energy margin is a convenient 
index for ranking the contingencies (3-phase faults at gene 
rator terminals). The proposed method can also serve as 

an accurate algorithm for contingency ranking. 



CHAPTER 5 


A STRUCTURE PRESERVING ENERGY FUNCTION INCORPORATING 
TRANSMISSION LINE RESISTANCES 

5.1 INTRODUCTION 

The structure preserving energy functions so far have . 
neglected the transmission line resistances. While this 
assumption can be justified due to high values of X/R ratio 
of extra high voltage transmission lines, it is used mainly to 
circumvent analytical difficulties of defining Lyapunov fun- 
ctions for lossy networks. 

With the reduced networks by assuming constant impedance 
load models, transfer conductances cannot be ignored and lot 
of effort went into the inclusion of transfer conductances into 
the energy function [30 , 50] . There is no adequate solution 
for this problem and the use of path dependent terms in the 
energy function is unavoidable. This can result in unreliable 
prediction of stability regions. 

The use of structure preserving model circumvents the pro- 
blem of transfer conductances and makes it possible to account 
for voltage dependent load models. Although arbitrary voltage 
dependent active power load can result in path dependent terms, 
it has been shown in Chapters 2 and 3 that the use of approxima 
tions by neglecting path dependent terms in the energy function 
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gives very accurate results. In general, it can be argued 
that path dependent terms which are relatively small in 
magnitude, can be easily handled by numerical integration. 

This makes it worthwhile, to investigate more detailed models 
of the system. 

A novel approach to develop a structure preserving energy 
function including transmission line resistances is presented 
in this chapter. This development of the energy function is 
based on the result that a network having the same conductance 
to susceptance ratio (G/b) for all its elements, can be trans- 
formed to a lossless network with a new set of power inje- 
ctions [67]. Even though G/B ratio in practical power system 
networks is not constant, it is possible to introduce approxi- 
mations by taking the average value of G/B. 

The structure preserving energy function based on the 
approximation is presented alongwith applications to the 
direct stability evaluation of two sample systems. 

5.2 NETWORK TRANSFORMATION 

It is shown in [67] that a lossy network with constant 
G/b ratio (g/B = a) for all its elements, is equivalent to a 
lossless network with a new set of power injections. The 
relationship between the two sets of injections is given by 
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Pi=P[ + o:C| (5.1) 

Ql = -a (5.2) 


where and refer to injections for the original network, 

Pj^ and Q]_ refer to the injections for the transformed (lossless) 
network. The susceptance matrices for both networks are identi- 
cal* Fig. 5.1 shows the original and the transformed networks 
with 'm* generators (internal) and 'n‘ load buses (including 
generator terminal buses) . 


From Eqs. (5.1) and (5.2), we can write 


p. = 1 . p - 

^i 7715 ^i 


1+a 


1+a' 


5 % 


(5.3) 


Q1 = 


a 


1+a 


5 Pi + 


1+a' 


5 «i 


( 5 . 4 ) 


The relationships given in Eqs. (5.1) - (5.4) are derived in 
APPENDIX I . 


5.3 -SYSTEM MODEL 

5.3.1 Generator Equations 

Consider the system shown in Fig. 5.1(b). Assuming classi 
cal model for generators, the motion of the ith machine can be 
written as (with respect to an arbitrary synchronous reference 
frame) 
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6 . = 

1 1 




6i = 0)- 


where 


mi mi 



(5.5a) 

mi'^ei* ^ ~ ^>2, m 

(5.5b) 

Eqs. (5.5) Can be written as 

(5.6a) 



il-Pgi > i = 1,2, ..., m 

(5.6b) 

“ Qei f i = 1,2, m 

(5.7) 


Eqs. (5.6) can also be written as (with respect to CX)I 
reference) 


where 


^i “ 


(5.8a) 




^mi ■" ^COI’ ^ •l,2,...,m (5.8b) 

^ "" l,2,...,m (5.9) 

bi[E|-E^Vi cos(e^-<P^)3,i=l>2,...,m (5.10) 


^COI = 


(5.11) 


In the above expressions^ bj^ is the susceptance of machine i 
The dynamics of the COI is described by 


^1^0 “ ^COI 


(5.12) 
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1 

where = jjq; _ M.u). 


m 

My = E M. 
i=l ^ 


5.3.2 Load Model 


(5.13) 


The real and the reactive power loads at bus i of the 
transformed network are modelled as arbitrary functions of the 
bus voltage, i.e.. 


^Li = i n 


(5.14) 

(5.15) 


5.3.3 Power Flow Equations 

The real arid reactive power injections at bus i of the 
transformed network can be written as 

n 

b.E.V. sin(<J),-e. ) + S B. .V.V. sin i = 1,2, .. . 

X J. i i 1 X J IJ 


PI = 


ni 


n 


2 ®xj^i^j ^-i' n » ^ ~ m+l, m+2, n (5.16) 


3=1 


^13 


cos(<Pi-e.)] - l ByV^Vj cos 

«3 


X 1. f2i f • 9 0 f in 


n 


2 B. .V..V. cos <^ 4 .:, i = m+l, m+2, ...,n (5.17) 
j=l ^ J 


where [b] is the susceptance matrix of the network 
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Since the net powers at bus i is zero, we have the 
following power flow equations : 


fpl = °> i = 1.2, ...,n 


(5.18) 


^i ^ ^qi ^^i^ i-i»2, ..,,n 


(5.19) 


5.4 STRUCTURE PRESERVING ENERGY' FUNCTION (SPEF) 


The system model described by Eqs. (5.8) - (5.19) is 

similar to the model considered in [54]. The only difference 

is that the transformed mechanical power P* . is not constant 

mi 

and varies with time. However, a general structure preserving 
energy function for lossless networks with time varying mecha- 
nical powers input is derived in APPENDIX C. ^^^plying this 
energy function for the system shown in Fig. 5.1(b), we get 


V/ = 


m 


1 2 
where = f 2 

i=l 


'•21= 


m t de- 
1=1 t„ 


(5.20) 

(5.21) 

(5.22) 



n 

Z 

i=l 


t d<}!. 

{ ^Lidr"^ 

^0 



n 

Z 

i=l 


Vi 

/ 

Vio 


Vi 


dV, 


(5.23) 


(5.24) 
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«i4 


lU p Q 

[E^+V^-2Ej^V^ «sO^^-(P^^))] 

(5.25) 


n n 




W^5 


5 ®ij( Vj ''ij-VioVjo ‘I’ijo) (5-2®) 


In the above expressions, subscript 'o* indicates the quanti- 
ties evaluated at the initial stable equilibrium point. 

Utilising the relations in Eqs. (5.3), (5.4), (5.7), (5.10) 
and (5.15), the energy terms Wil- W ^2 ^23 written as 

m t 2 de,- 

*21 = ■ 1=1 { wr 


m t d0. m t « 

■ { I’mi 3t" <11 + cos(Sj-<Pj) 


d9, 


d9^ 

dt 


dt 


m 
Z P 
i=l 


m t „ de. 

mi^ V\o> “ { "i^^i-Vi <^°3(9,-<<P,) 


dt 


- + W2g 


m 


where, «^i = - . P„j. ( V®io> ' 


m 


t ^0. 

™26 “ dT ®1 

i=^l 


(5.27) 

(5.28) 

(5.29) 
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W 




22 


n t , d<p. 

2 S t““ 6 rTF” 

i=l t 1+a 
o 


(5.30) 


V. 

n 1 


W, 


23 


= 2 / (• 
i=lV.^ l+a 


aP. .+ 




V^)3 dV. 


(5.31) 


The energy function W can now be expressed as 


W = = w, . w,, 


(5.32) 


where 


^22 “ ^22 * ^23 


^23 * ^ 24 - ~ ^24 ^25 ~ ^25 * 


Comments 

1) As described in Chapter 3, if the mode of instability is 
known apriori, the kinetic energy term can be replaced by 
Wii which corresponds to the kinetic energy coirponent res- 
ponsible for system separation. If it is assumed that the 
faulted generator k separates from the rest of the generators 
(which swing together) » then 


W 


11 


1 ^T-k 


f“k ■ “T-k>' 


(5.33) 
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2) The energy term in Eq. (5.29) can be written as 


'li. 


26 


=2 f 


^ dQ, 


r ifCp i ^ 

= [eb.E^ / g_ dt] - [cb.E. / cos{9^-<p^) ^ dt] 

^2 ni t rin ■ 

= ab^E^(0^-0^^) - I [ab^E. / V. cos(e.-<P^) ^ dt] 


^261 ^262 


v/herc % 


m 


'261 = “Vi^°i-®lo> 


(5.34) 


(5.35) 


W, 


262 


m t de. 

ilj_ { ''i 3r ‘5-36) 


In tine above expressions, while the term W 2 ^j_ is path inde- 
pendent, the term ^^^2 dependent and can be evaluated 

by numerical integration using trapezoidal rule [75]. 


3) For constant real and reactive power loads, the term W, 


andWl 23 can be directly computed from 

“ ih fPLi-“«Ll>‘''i- ''io>3 




22 


(5.37) 


(5.38) 


1 1 is interesting to observe that an additional condition of 
reactive power components of the loads being independent of the 
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2) The energy term (5.29) can be written as 

m 2 /. t . dg 

^26 ^ .L [ah^Ei f ^ dt] - E [ab^E^ / cos{ 0 .-<p^) ^ dt] 

X i 1 X 

0 

m 2 m t do. 

~ "* cos( 8 ^-<P^) dt] 

X "“1 1*~1 t 


- + ^^262 


(5.34) 


m 2 

where 


(5.35) 


W, 


262 


m t dS. 

S ab.E. / V. cos( 0 .-<P.) dt 

i_l ^ ^ t ^ 1 1 ' dt 


(5.36) 


In the above expressions, while the term ^ 2 ^^^ path inde- 
pendent, the term '^^ 2(^2 Ps^'^h dependent and can be evaluated 
by numerical integration using trapezoidal rule [75]. 


3) For constant real and reactive power loads, 
and W 23 can be directly computed from 


*22 = 


“aa = V^°9e 


the term W 22 


(5.37) 


(5.38) 


It is interesting to observe that an additional condition of 
reactive power components of the loads being independent of the 
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voltage, is required in this case compared to the case in loss- 
less network to make the terms ^22 ^^23 independent. 

4 ) The component W 2 £, is always path dependent as 


m t d©. 

^262 = { ''i ar 


(5.39) 


cannot be integrated exactly. This implies that there is 
always a path dependent term in the energy function when the 
transmission line losses are considered. However, if this 
term is small in magnitude, it should not affect the accuracy 
in prediction significantly. 


5) For practical power systems, G/B ratio for all the elements 
is not the same. However, the proposed energy function can 
still be employed to predict the critical clearing time using 
an average value for G/B. 

PEBS method is employed again, as in Chapter 3 for a loss- 
less system, to evaluate the critical energy for a fault at the 
terminals of a generator in a lossy network. The computer 
program described in Chapter 2 can be used with a minor modifi 
cation to incorporate the calculation of the energy term ^26 
in Eq. (5.29) in the computation of SPEF. This term is path 
dependent and can be evaluated by numerical integration. To 
accommodate the known mode of instability, the kinetic energy 


term is also modified. 
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5.5 numerical examples 

5.5.1 Description 

Two system examples (7-machine and 10-machine) are 
studied to investigate the transient stability of lossy power' 
system networks. The single-line diagrams alongwith system 
and operating data are given in .APPENDICES E and F. 

Three phase faults at different buses are considered for 
analysis. Fault clearing is followed by instantaneous reclosure 
of the line. This implies that the prefault and the post- 
fault network configurations are the same. Two different load 
characteristics are selected to illustrate the applicability 
of the proposed method to voltage dependent load models. 
Different load buses are assumed to have identical characteri- 
stics for the sake of simplicity in analysis. The following 
cases are investigated. 

Example 1 : 7-iiiachine 10-bus (CIGRE) test system 

a) Fault at bus 1 ; b) Fault at bus 2 } c) Fault at bus 3; 
d) Fault at bus 5 j e) Fault at bus 6 ? f) Fault at bus 7. 

Types of load characteristics t 

i) kj^V^ + jk^V^ J ii) 


Average a =0.2. 
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Example 2 : 10-machine 39-bus (WEW ENGLAND) tek system 

a) Fault at bus 1 ; b) Fault at bus 3 ; c) Fault at bus 4 ? 

d) Fault at bus 5 j e) Fault at bus 6 ; f) Fault at bus 7 ; 

g) Fault at bus 8 ; h) Fault at bus 9 . 

Types of load characteristics : 

i) + jk^V^ ; ii) 

Average a = 0.07. 

In both the cases, results are also obtained for two 
arbitrary values of a } (i) 0.1 and (ii) 0.2. 

5.5.2 Results 

In all the cases, critical clearing time is predicted 
using the proposed energy function for the assumed and average 
values for G/B and load characteristics of types (i) and (ii) . 
The result is compared with that obtained by digital simulation 
and presented in Tables 5.1 - 5.2 and 5.4 - 5.6 for 7-machine 
and 10-machine systems respectively alongwith the critical 
energy. The predicted results are quite close to the results 
obtained by digital simulation. 

Critical clearing time is also obtained by digital simula- 
tion for the average value for G/B of the systems under investi- 
gation and their actual line resistances for load representa- 
tion of types (i) and (ii) and compared in Tables 5.3 and 5*7 
for 7-raachine and 10-machine systems respectively. 



Table 5.1 Effect of transmission line resistance on critical energy,'// and 

' cr 

predicted critical clearing time, t._ for 7-machine System 
a = 0.1 ; in sec and W^p in pu 


178 


CO 

> O 

CM +> 

•f— ) I 


m o o in :)c 

CM CO O' 

'sT CO CO I 

• • • « 

oo oo 


O UT) He 
OO 
vOvO I 

' * » 

OO 



vO 

* 

5i< 



00 



o 


00 

1 

1 

\0 

* 

• 



» 

o 

o 



o 

1 

' 

1 ^ . 

o 


* 

iD 

I ^ 

CO 



o 

. t;:}- 

CO 

1 

1 

iC 

• 

• 



• 

O 

o 



o 


(0 1> 
O vO 

•H 0> lO 

^ u r- 

•H O U • 
C O 00 
O C14B: rH 


> M 

H U 

M -P 


T5 iH 
03 03 


C • 

O O 
O 2 Jd 


I I I I I I 

Oif) oin Oio in o Oin Oin 

CT^ C3^ uO iT) vO \D HI 04 H *H vD \0 

00 00 Kj- ^ vOvO COCO 

• • •* ••'•4 •« 

OO oo Oo Oo oo oo 


in 

\0 

H 


CD 

H 

00 


vO 

H 

O 

\0 

00 




vO 

00 

* 

• 

• 

• 

• 


o 

o 

o 

o 

o 

o 


o S*' 

•H D» 

■p (-f 
•H O U 

h 

U PU 3: 


00 

CO 

0- 

00 

r- 

O 

00 


lO 

H 

o 

vO 

CO 




vO 

00 

• 

• 

* 

• 

• 

• 

o 

o 

o 

o 

o 

o 

lO 


vO 

CO 

lO 

l> 

CM 

0- 

H 

LO 

00 

CM 

lO 


CM 


CM 

vD 

CO 

H 

H 

lO 


lO 

« 

• 

• 

• 

• 

» 

in 

vO 

0- 

o 

00 

MO 

iH 

iH 

H 


H 


rH 

CM 

CO 

lO 

vO 

r- 






C 



o 


c 

*H 


o 

P 


•H 

O 


-P 

•H 


O 

•o 


•H 

CD 


TO 

P 


CD 

a 


M 



a 




Si 


s 



n 




p 


>• 

•H 


-p 

H 


•rl 

•H 


iH 

JO 


*H 

CO 


JO 

p 


CO 

CO 


■P 

c 


if) 

•H 


C 



♦H 

P 



O 





O 


c 


0 

o 

Q) 

T* 

•p 

TO 

O 

p 

O 

S 

CO 

e 


H 


0 

D 

o 

x: 

a 

x: 

p 

•p 

-p 


(0 


CO 


C3> 

c 

rP 

c 

•p 

<0 

•H 

p 

p 

-P 

0 

•p 

O 

T3 


o 

•rl 

•P 

H 

<0 

TO 

D> 

C 


G> 

O 


C 

o 

X> 

T5 

TO 

T5 

O 

0 

0 

C 

C 

C 

•H 

•P 

•P 

CO 

0 

0 

+> 

P 

P 

iO 

JQ 

JO 

o 

O 

o 

CD 

0 

0 

c 

J3 

a 

•H 

P 

•p 


P 

p 



D> 

C 

C 

C 

•H 

^P 

•P 

P 

P 

P 

CO 

CO 

«t3 

O 

0 

0 

rH 

rp 

fP 

O 

O 

O 

H 

rP 

H 

CO 

CO 

0 

O 

o 

O 

•H 

•p 

•P 


p 

P 

•H 

•p 

P 

P 

p 

P 

o 

O 

O 

H 

II 

II 


CM 

CO 

P 

P 

p 

U 

O 

u 

P 

P 

p 








in 

o 



lO 






00 

Ph 


0 

p 

♦ 



CM 

vO 






o 

p 


0 

1 

o 

in 

* 


CD 

{ 

o 









o 

1 

1 

CM 

1 







• 

'M' 

• 


\0 







o 

o 



6 





CM 



c> 

p 



rr^ 





•r-> 

o 

p 


CO 

« 

o * 

• t 

* 

1 

CM Jfc 
vO 

• 1 

X5 

C 



p 



o 

o 


o 

1 

CO 



rH 












-id 









u 

( 



p 

1 

t p 


f 






o 




1 m 




00 





•H 


CO 

O' sfe 

00 

* 

xj ^ 

vH 


S 


• 5 


* 

« 

1 

1 

• 

1 

>- 

cn 

CD 

P 


• * 

fH 



; o 

o 

o 

M 

CO 


CD 



1 







>. 


T5 1 



1 






c 

(0 


O 









<D 



S 

p 









<D 


XJ 

CO 


sO 

o 





iH 

c 

1 

o >• 


O 

CO 



00 


<0 

O 

•HI 

x: 


(0 

o 

•H Cn 

p u 


lO 

(M ♦ 

00 

5ic 

tf) 

ir> 

*rl 

o 


p 

P CD 

u 

« 

« 

1 

1 


1 

-P 

(0 



u c 

o 

vO 

vO 

NO 

•rl 

tj 

S 

1 

C-- 



OWB: 

p 

p 



P 


M 

u 



1 








C 

u 





o 

in 

in 

o 

o 

lO 

o 

o 



I 


p 


CO 

CO 


00 

<D 

p 



CO 



« 

0 

• 

vO 

• 

CO 

• 

O 

M 


<n! 

1 M 


o 

o 

o 

o 

o 

o 

C 

O 

ii 

o 


1 

1 

1 

I 

1 

CO 

p 


•o 1 

p 


lO 

o 

o 

in 

in 

o 



Z3 

J 



o 


00 

CM 

CO 

00 

(0 

«v 

CL 

+ - 



'M- 


'M- 


nO 

CO 


o 






0 

* 

'0 

0 

• 

<0 

0) 

S 

in 

^ . 



o 

o 

O 

o 

o 

o 


P 


P 









0) 

D> 

u 

M 

CM 


CO 

1> 

CM 

in 

CO 

o 

o 


u 


O 

vO 

00 

CM 

CO 

t- 

c 

C 



o 






vO 

CO 


•H 


•H 

p 


0 

• 

0 

• 



H 


*o 




o 

o 

o 

o 

o 

o 


CO 

C 

0* 









C 

<D 

CO 










O 

P 


P 









•H 

O 

o 

0> 

p 


p 

in 

o 


p 

00 

0) 


o 

X5 

u 


Q 

vO 

00 

CM 

CO 

r- 

W 

rH 

CO 

O 

o 






\D 

CO 

•H 

(0 


s 

p 



* 

0 

♦ 

0 

* 

6 

o 

c 




o 

o 

o 

O 

o 

o 

<0 

•H 

•H 

X3 









c 

P 


CO 









<0 

•H 


O 

p 









U 

o 

P 

CO 



o 

CM 

00 

vD 


•p 

O 

p 


o 



r- 

vO 

00 

O 

o 





•H CO 


p 

so 

P 

(O 

P 

r- 

X? 



P M 


00 

sO 


vO 


in 

o 

0) 



P 0 


• 

• 



0 

• 


P 

CM 


M c 

o 


CO 

o 

o 

in 

r- 

p 

o 

# 


O W3: 

p 

p 

p 

p 

p 


o 

•H 

O 










0) 

X> 











p 

(D 

II ■ 

' •' 4 

► 








p 



C O 


p 

CM 

CO 

in 

vO 

c^ 

w 

a 

55 

<u z 










CD 









<N 












« 



c 









lO 



0 

> 











•H 

1 








0) 



PP 


p 

CM 

CO 

in 

nO 

r^ 

rH 



P CO 








JQ 



CJ o 


0> 

U) 

if> 

CO 

<0 

0) 

CO 



<0 o 


D 

O 

o 

o 


o 

H 





CQ 

CQ 

DO 

CO 

CQ 

CO 


I 





179 



Table 5.3 Comparison of critical clearing time, t by digital simulation for 
7-machine system 
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Table 5,4 Effect of transmission line resistance on critical energy, and 

predicted critical clearing time, t^^, for 10-machine system 
a = 0,1 ; t^p in sec and W^p in pu 
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Table 5,6 Effect of transmission line resistance on critical energy, VV and predicted 

C JT 

critical clearing time, t^^ for average value of a for 10-machine system 
Average value of a = 0.07 ; t^.^. in sec and in pu 
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Swing curves for examples 2-a)(i), 2-b)(i), 2-c)(i) ,2-f ){i) 
and 2-g)(i) are shown in Figs. 5.2-5.11 for stable and 

unstable cases for G/B = 0.1 and load characteristic of type (i) . 

The variation of total energy, kinetic and potential 
energy components for examples 2-a){i), 2-b) (i) ,2-c) (i), » 2-f) (i) 
and 2-g)(i)are shown in Figs. 5,12-5.13, 5.16-5.17, 5.20-5.21, 
5.24-5.25 and 5.28-5.29 for stable and unstable conditions for 
the assumed value for G/B = 0.1 with the loads modelled as 
type (i) . The corresponding variation of the components of 
potential energy function, W2 are shown in Figs. 5.14-5.15, 
5.18-5*19, 5.22-5.23, 5.26-5.27 and 5.30-5,31 respectively. 

The variation of the components of the energy term is shown 
in Fig, 5*32 for the stable case with G/B - 0.1 and loads 
modelled as type (i). 

5.5.3 Discussion 

It can be observed from the tables that the proposed 
energy function including the transmission line resistances 
yields a higher estimate of the critical clearing time for 
higher values of resistances. Also, the predicted values are 
in close agreement with the results obtained by digital simula- 
tion. 

4 

From Tables 5.3 and 5.7 it can be seen that the use 
of an average value for G/b gives accurate results in con^arison 



Rotor angles ( rad) > Rotor angles (.rad) 
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FIG. 5,2 SWING CURVES FOR 10-MACHINE SYSTEM, STABLE 
CASE 


Roto 



FIG.5.3 swing curves FOR 10-MACHlNE SYSTEM^ 
UNSTABLE CASE 



angles (rad) 


-Foult cleared at 
0.240 sec 


Fault at bus 3 
Example 2-b)(i) 

Load Model:!) KiV^+j K2V^ 
0/B = 0.1 




Time m seconds 


-Fault cleared at 0.240 sec 



1.0 

Time in seconds 


FIG. 5.4 SWING CURVES FOR 10-MACHlNE SYSTEM, STA BLE 
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Fault at bus 3 
Example 2-b)(j) 



FIG. 5.5 SWING CURVES FOR 10- MACHINE SYSTEM, 
UNSTABLE CASE 



Rotor angles (rad) ^ , Rotor angles (rad) 


190 


Fault at bus 4 
Example 2-c ) (i) 




FIG. 5.6 SWING CURVES FOR lO-MACHlNE SYSTEM, STABLE 
CASE 
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PIG.5.8 SWING CURVES FOR 10-MACHINE SYSTEM, STABLE 
CASE 
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FIG.5.9 swing curves FOR 10-MACHINE SYSTEM, UNSTABLE 
CASE 
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STABLE CASE 



Rotor angles (rad) Rotor angles (rad) 
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FI6.5.13 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 
FOR 10-MACHINE SYSTEM, UNSTABLE CASE 
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FIG. 5.14 VARIATION OF THE COMPONENTS OF WoFOR 10 
MACHINE SYSTEM, STABLE CASE 


^Fautt cleared at 0.275 sec 



Time in seconds 


-40.0‘ 


FIG.5.15 variation OF THE COMPONENTS OF W 2 FOR 10 
MACHINE SYSTEM, UNSTABLE CASE 
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Fault at bus 3 
Example : 2-b) (i) 
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Fault at bus 3 
Example : 2- 

Load Model : 
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/" W21 


.40oLp>G.5.18 variation OF THE COMPONENTS OF ^2 FOR 
10-MACHlNE SYSTEM, STABLE CASE 

40.0r 


20.0 


-Fault cleared at0.245sec 
rW24*VV25 
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FIG.5.19 VARIATION OF THE COMPONENTS OF W 2 FOR 
10-MACHINE SYSTEM, UNSTABLE CASE 
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Example ; 2-c) ( i) 
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FIG. 5.22 VARIATION OF THE COMPONENTS OF W 2 FOR 
10 -MACHINE SYSTEM, STABLE CASE 
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FIG. 5.23 VARIATION OF THE COMPONENTS OF W 2 FOR 
10 MACHINE SYSTEM^ UNSTABLE CASE 






£,\JL 


Fault at bus 7 
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FIG.5.26 VARIATION OF THE COMPONENTS OF W 2 
FOR 10- MACHINE SYSTEM, STABLE CASE 


50.0 



F1G.5.27VARIAT10N OF THE COMPONENTS OF W 2 FOR 
10- MACHINE SYSTEM. UNSTABLE CASE 
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FIG. 5.29 VARIATION OF TOTAL ENERGY AND ITS COMPONENTS 
FOR 10-MACHINE SYSTEM, UNSTABLE CASE 
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FIG. 5.31 VARIATION OF THE COMPONENTS OF W2 FOR 
10 -MACHINE SYSTEM, UNSTABLE CASE 


W26 
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i 

i Foult at bus 1 


Example 2-a)(i) 

Load Model ; i) KiV^* j K2V2 



FIG. 5.32 VARIATION OF THE COMPONENTS OF W26FOR 
10-MACHlNE SYSTEM, STABLE CASE 
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to those obtained considering the actual line resistances. 

Thus* the proposed method to account for the resistances is 
quite satisfactory. 

The replacement of the kinetic energy term by 
results in a slight improvement in the predicted result of 
critical clearing time. As in Chapter 3, was calculated 
by considering the faulted generator going out-of-step with 
respect to the rest of the generators. As observed in 
Chapter 3, this procedure of the kinetic energy correction 
does not significantly affect the results. 

Out of the 6 components of the potential energy, ^^ 21 * '^24 
and W25 are identical to the similar components defined in 
Chapters 2 and 3 and their variations are similar to those 
observed in Chapters 2 and 3. It is also observed here that 
w,, is significantly larger than 1*23. However, in comparison 
to’ the lossless case, the magnitude of W23 is slightly larger. 
It is interesting to observe that in all stable cases W2g^ 

(which is positive) is of similar magnitude as of ^23 (which is 
negative) and tends to cancel each other. In unstable cases, 
^26 increase without limit although the rate of 

increase is smaller than the rate of decrease of 

The variations of the component ^ 2 ^ alongwith its sub- 
components W25 j^ and W252 shown in Fig. 5.32 for example 
2-a)(i). It is interesting to observe that the path dependent 





term W252 negligible initially and also which is path 
independent, dominates over This shows that the accuracy 

in the prediction of critical clearing time is not significantly 
affected by the presence of the path dependent term whose con- 
tribution is negligible. 

5.6 CONCLUSIONS 

In this chapter, a structure preserving energy function 
has been developed including the effects of transmission line 
resistances. The function has been applied to two sample power 
systems and based on this study, the following conclusions are 
drawn : 

i) The predicted value of t^^ agrees closely with that 
obtained by digital simulation. 

ii) Inclusion of transmission line resistances in the energy 
function gives an estimate of t which is, in general, 
higher than that obtained for the lossless system. 

iii) t obtained by digital simulation using the average 
value for G/B agrees well with that obtained using the 
actual resistances of the elements of the system. Thus, 
the approximation introduced by the assumption of same 
G/B ratio for all the elements of a network, appears to 
be fairly accurate. 



CHAPTER 6 


A STRUCTURE PRESERVING ENERGY FUI'ICTION WITH DETAILED 

GENERATOR MODELS 


6.1 INTRODUCTION 

A structure preserving energy function was developed in 
Chapter 3 representing the generator by the simple model, i.e. 
a constant voltage source behind the transient reactance. 
Recently, Padiyar et al [58] have considered detailed models 
for generators including flux decay, transient saliency, AVR 
and exciter to evaluate the transient stability of power 
systems. The SPEF considered in [58] has a time derivative 
that, along the post-fault trajectory, is non-positive which 
results in pessimistic estimates of the region of stability. 

The principal aim of this chapter is to reexamine the 
work reported in ref. [58] and try to remove the limitations 
mentioned above. A new structure preserving energy function 
is developed which is constant on the post-fault trajectory. 
This should lead to better accuracy in the prediction of 
stability using PEBS method. Moreover, the generator model is 
improved upon by the addition of a damper winding on the qua- 
drature axis. Simpler and more general expressions for the 
energy functions are then derived where the use of the classi- 
cal model becomes a special case. The simplification in the 
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expressions can lead to increasing the speed of computation 
of SPEF. Investigations are carried out on a 4-machine and a 
lO-machine system examples to illustrate the proposed method. 


6.2 SYSTEiVl MODEL 


6.2.1 Generator Model 


Consider the 'm’ machine system supplying 'n' nonlinear 
voltage dependent loads, as shown in Fig. 3.1 (Chapter 3). 
Usually, in transient stability studies, the generators close 
to the fault are modelled in greater detail [24, 57]. The 
complexity in modelling the generators more realistically 
depends on the number of rotor circuits represented. The 
model considered here includes two circuits on the rotor; 
i) the field winding on the direct axis and ii) a damper 
winding on the quadrature axis. It is assumed that the mecha- 
nical input to the generator is constant. The dynamics of the 
machines are then described by the following differential- 
algebraic equations (with respect to 001 reference frame) ; 


For i = 1,2 , . . . , m 





M 


i 


CO. 

1 




^COI 


( 6 . 1 ) 

( 6 . 2 ) 



where 

Si * ®dl idl 
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(6.3) 

m 



^COI "" ^^mi~^ei^ 


(6.4) 

^doi ^qi = ^fdi ~ ^qi + 

•Si> Idi 

(6.5) 

^qoi ^di " “^^di ■*“ ^^qi"^qi^ 

Si^ 

(6.6) 

^qi ^qi ^di ^di 


(6.7) 

^di " ^di “ ^qi ^qi 


(6.8) 

e "= (V^, + jV,,) e ^ 


(6.9) 


From Eq. (6.9) we get, 


^qi “ 


( 6 . 10 ) 


^di ~ ‘■'^i sin(e^ - <P^) 


( 6 . 11 ) 


From Eqs. (6.7) - (6.11), and are obtained as 

^di =[V. cos (e. -dj) - Ey/x;j. (6.12) 

(6.13) 


^ _ E^i + slnCe^-d^) 


qi 


qi 


Substituting for and i^^ from Eqs. (6.12) and (6.13) in 
Eq . (6.3 ) and simplifying, we get. 


^ei 


“^oi ^i sin(ei-<P.) E^j. V. cos(0i-<Pi) 
— ■ + 


= -ai 


X 


di 


qi 


+ V? sin 2 (0.^.) (6.14) 

^ 2xJjiXj. 

Eliminating i^^ and i^^i^ innqs. (6.5) and (6.6) using Eqs. 
(6.12) and (6.13), we get, 



T* F' 
^doi ^qi 


T' F' 
^qoi ^di 


X 

E. .. - -Si F* + 
idl y* “«■? ^ 
Xdi qi 


Vj cos(Q^-qn,) 


Y • 

^di 


**di~*di* 


= -[> E. . 


•X’. “di 
qi 


qi 


(x ^-x' . ) 
qi qi^ 


(6.15) 

(6.16) 


6.2.2 Excitation System and Voltage Regulator Model 

The excitation system and voltage regulator model used here 
is IEEE type 1 taken from ref. [76]. The block diagram of the 
system is shown in Fig. 6.1. The system is described by the 
following dynamic and algebraic equations : 


^ei 

~ ^Rfifi ■* ^i - '^si 

(6.17) 

« 

''ri 

= ''ei - V)Aa 

(6.18) 

VRLi 

VRi if i Vr^ < 

(6.19) 


’^Rmin ^Ri ^ ^Rmin 

(6.20) 


~ ^Rmax \i ^ ^Rmax 

(6.21) 




^fdi 

r 

= tVRLi " A,, e - Kg E,^,]/Tg 

(6.22) 


K 


'^si 

“ tJ ^fdi " ''^sLi 

(6,23) 


where, 


^sLi 


K 

(^ E 


T 

s s 


fdi “ ^sLi^ 


(6.24) 
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6.2.3 Load Model 

The load model given in Chapter 3 is used. 

6.2.4 Power Flow Equations 

For the lossless system shown in Fig. 3.1 (Chapter 3), 
the following equations can be written at bus i : 


Pi = 


Egl Ejil \ OOs(et-»t) 


Y t 

^di 


X' . 


Vf (x’ -X’ ) 

+ -i — — c[i_ 2(<p.-e.) 

1 1'^ 




Qi 


2x 

1 

di 

x' . 
qi 

n 

£ 


V,- 

j=l 


1 

n 

E 

B. . 

V. 

3=1 

13 

01 

1 

''i 


•f • • • > 


m 


E B. . V. V. sin ^ , i = m+l,m+2, ...,n (6.25) 

j =1 ^ J 

pt V eiri/’/n _a \ 
c^i si 


1 1' 


1 1' 


^di 


X 


di 


’'qi 


n 


• “"X * • ) 

i — di — c[i_ [(,q 32 (0.-(p. )-1] - Z B. .V.V. cos <P. . , 

Oyl v' ^ ^ 

^^di ^qi ^ 

1 ~ m 

n 

= - Z B. . V- V. cos i = m-i-l,m+2, •• • ,n {6.26} 

^ J 

In the above expressions = Im where [Y] is the 

admittance matrix of the network (excluding machine reactances) 
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The power flow equations, as given in Chaliter 3, are 


Pi fpi(V.) , 


Qi-VV • 


1 ~ 1,2, •••, n 
i ~ 1,2, n 


(6.27) 

(6.28) 


6.3 STRUCTURE PRESERVING ENERGY FUNCTION (SPEF) 

In what follows, a structure preserving energy function, 
which incorporates the effect of AVR, transient saliency, damper 
winding and voltage dependent active and reactive power loads, 
is developed. 

Consider the function defined for the post-fault system, 

13 

W(e,u),E' ,EL,V, <P ,Ef^,t) = Wj^+W2 = Wj^ + 2 W2i (6.29) 

4 1—1 


where 


W^(u)) = 

1 ^ 2 

4 2 M.ojf 

^ i=l ^ ^ 

(6.30) 

W2i(e)= 

m 

- 2 P„4(0:i-9iO^ 

^ mi 1 10 

(6.31) 

II 

n t dtp 

2 / ^pi^^i^ dT 

i=l 

(6.32) 

W23(V)= 

n 

E / dV, 

i=l ''i 

(6.33) 

W24(®'^ 

' = ? [EL? + V? -2E'iV. cos(9-<Pi) 

q 1=1 4 



(6.34) 
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ih »ijo) (6-35) 

m 

W 24 {e.V,<P)= - (cos 2(e^-p^)-l) -v2_^(cos2(e. -„. )-l }] 


m t 

=-z / 


X ’ j .- X ’ . 
01 ai 

4xL. x» . 
01 qi 


‘fdi 


■‘^0 


dE' . 


m E’? - E'? 

W . o(EM - E -ai ai 2 

a 2( X -x' ) 


W 29 ( t ) 


m 


t TL., 


Z / -: doi (^ 2 


'o '"di ^di^ 
m 


'2io(e.Ei.V,<f) =^£^ [Eif+vf+2EJlVj sin(9.-^j) 

- t ) 1 -^ 


qi 


’211 


’212 


(V) 

m 

= - E 


i=l 

(t) 

m 

= E 


i=l 

(Ei) 

m 

= E 

• —I 


V? - V.^ 
1 10 


'qi 


T' . 
qoi 


< 5 t ^)" dt 


o ' qi qi- 

: t 2 gt 2 

'di ^dio ■ 


(6.36) 

(6.37) 

(6.38) 

(6.39) 


(6.40) 


(6.41) 


(6.42) 


(6.43) 


Jt. — u. V ^ ^ Z J 

qi qi 

Subscript • o’ in the above expressions indicates the quantities 
at the initial sep. 



Now it will be shown that the time derivative of the 
function defined in Eq. (6.29) is zero. 

Partially differentiating % with respect to E' . , E' . , 

^ qi di ’ 


dW^ 

2 

qi 


E’ . 

-Si . 

y t 

^di 


V. cos(e.--9).) E'. 

i — i_ ^ — ai. 


X’ 

dr 




'='di-=‘di) 

E* . 


x],. 

di 


E 


X 


(’‘dr^di) 


- _ doi qi 


fdi 


^^di“^di^ ^^di"^di^ 


, from Eq. (6.15) 


(6.44) 


dV(2 

A'p I 

^‘^di 


Edi , ^1 

(Xql-qi) 


= ial 

X’ 


E^l ^ V. sln(6.-d^) 




qi 


- laSiik , from Eq. (6.16) 


(x . -x’ . ) 
qi qi 


W7 


(6.45) 


, V? E'.V. cos(e.-V.) n 

-L.\f (v )+ i— - Z B. .V.V. cos <P. . 

V. •-^qi'' x^^ x^. ^=1 ij 1 3 13 


^"di^ii 


i { cos2(0^-4).)-l} + 


di 3=1 

E',V, sin(6^-dt), 


^ f i— 1 ,2 » . • • »ni 


qi 


= ^ ff .(V. ) - Z B. dViV. cos J> i = m+l,m+2,..., n > 
V. •- qi' 1 '' •_•, 13 1 3 13 

X J 


n 
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1 

V 


- i:Qi+fqi{Vi)] , i = 1,2, ..., „ 

in view of Eq, (6,26) 


= 0 , from Eq. (6.28) 


(6 .46) 


dW2 

■PT 


dV(2 

^0. 


^qi \ sin(0^~(p^) n 

V?(x’ -X’ .) 

- 3^ di ' 


oi^ . . V cos(0.-<p,) 

— sin2C0.-^.) 1 ^ ^ i=:l 2 

, 1 l' x' . 


2x' . X* . 
di qi 

n 

E B. .V.V. sin (p. . 
j=l ij 1 3 ij 


qi 


1 ■“ m“i“ 1 f m^"2 y y o 


= » from Eq. (6.25) (6.47) 

sin(e.-<P.) V^(x^i-x-i) _ 

+ — + a — 53 _j 


‘di 


^^di ^qi 


sin 2(0^- <P^) 


E^i V. cos(9,-.,) 

qi 


= + Pgj^>in view of Eq. (6.14) 


(6.48) 


dW. 


n 


d<P^ 


m E 




fdi 


i=l 


m 
+ E 


T» 

Moi 


•=i (>‘di-’‘di) 

T' 
goi 


Si 

dt 


1=1 


dE' . m 

( j+r^) + 2 

i=l (x 4 -x'.) 
' qi qi 


"Edi.2 


(Ht^)" (6.49) 


Partially differentiating Wj^ with respect to we have. 


dWj^ 

— M* to . 
00 ). 1 1 
1 


.m 


(6.50) 



219 


Substituting from Eqs. (6.48) and (6.50) we get, 


dW2 de^ 

ss;; ar + 35;; dT = - Pmi“i + Pei“i 


= (M^-ol. - + p^.)^. 


“i^i 

Poor from Eq. (6.2) 


(6.51) 


Thus, 

m dW, do). dvIL d©. 


^3S7 ar * SeT ar = - 


m 


^COI 


^COI 

S M.w. = 0 , from the definition of CX)I 


M 


T i=l 


3. 1 


variables 


n dW 2 dV^ 


2^ dft” ~ ^ * from Eq. (6.46) 


(6.52) 

(6.53) 


Substituting from Eqs. (6.47) and (6.49), we have, 


n dW 2 d<P^ dW 2 n 


m 


i!i dr> ar = 


Exj- dE’ , 
fdi gi 


^^di"^di^ 


dt 


m TL . 

+ r 


i=l 


Mt ' * ,i, xMt ' 


i=l (x .-x' . ) 
' qi 


3t ^ J, v^dt > ■ 


= - z 

i=l (x,,-x' ) 


i-1 (Xjjj^-X^^) 


di ^'di- 

+ 2 — - (;t+^)^» E q.(6.27) 

(6.54) 




‘dt ' 
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Substituting from Eqs. (6.44), (6.45) and (6.54), we get, 
m dWo dE' . m dw« dE!,. n dW^ d<p. d% 

'*4 ^ 

; i. ‘ui *■! ] , ; . 

i=l ^^di“^di^ ^^di“’^di^ ^^qi”^qi^ 

gfdi ^1,1 ^ I 


= 0 


(6.55) 


Hence, 


dt “ dt 


{'ft 1+^2'^ 


m 

= 2 
i=l 


dWn do). SVIU dS. n 

ar + 557 ar) + 

1 1 i-J. 


dWo dV. 

(s^ ar5 





= 0 ,in view of Eqs. (6.52), (6.53) and (6.55) (6.56) 


Comments 

1) The last 4 terms W2io» ’^211' ^212 ^213 correspond to 

the change in energy due to the presence of damper winding on 
the quadrature axis. In the absence of this winding, E^ = 0 and 
these terms vanish. Thus, there are only 9 terms in the energy 
function when one-axis model of the generator is used. Out of 
these, the first 8 terms are the same as those given in [58, 71]. 
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The 9th term is introduced as the energy function defined here 
is constant along the post-fault trajectory. It is to be noted 
that the energy function defined in [58,' 71] has a time deri- 
vative given by 


- Z 


which is non-positive. 




(6.57) 


2) The term W 2 ^ in Eq. (6.36) is the change in energy stored 
due to saliency in absence of which it is zero. 

3) In the absence of AVR, is constant and in Eq. 
(6.37) can be expressed as 


_ 2 E . 

i_l fdi 


(6.58) 


With this simplification the sum of the first 8 terms of the 
energy function are equivalent to that given in ref. [57] for 
the one-axis model including flux decay but neglecting AVR. It 
is to be noted that ref. [57] considers only constant P-Q type 
loads. 

4) If a fault at its terminals separates generator k from the 
rest of the generators, then the kinetic energy term can be 
modified to W 2 ^j_ defined below (see Chapter 3) 




1 ^ k^T-k / \2 

M— ^‘^k-^T-k^ 


(6.59) 
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5) As discussed in Chapter 2, the potential energy compo' 

nent due to active power loads, can be approximated as- 

m 


- 2 [f . (V. ) <P. - f .(V. ) <P. ] 

pi 1 1 pi'' lo' ^io-* 


(6.60) 


A Simpler Expression for SPEF 

The terms ^ 24 * ^’’^25' ^26’ ^^210 

change in energy stored in the machine reactances and transmi- 
ssion lines. It is shown in Appendix J that this change is half 
the sum of the change in reactive power loss in the network, 

1 • 0 • > 


k ^OSS 

. m n 

= 5- [ ^2^ (Qai-Qaio) ’ ^^i~^io^^ 

From Eq. (6.5) we get, 


i 


di 


^fdi ^ ^gi ^ ^doi ^gi 

(Xdi-x*i) (Xdi-x'.) (x^i-x'.) 


(6.61) 


(6.62) 


Multiplying both sides of Eq. (6.62) by dE^^/dt and integrating, 
we get the following equation after summing over all values of i 


m t 

2 / i 

i=l t 

o 


dE' 
di dt 


dt = - 2 


m t 


■fdi 




dE', 

cit^ 


m ^ E' . dE’ . 

+ . \h — , T dt^ 

i-l^Xdj^-Xdi) 

^ ^ 32i_(^)2dt 

= Wo^+WoQ+WoQ » in view of Eqs.(6,.37 )-(6 .39 ) 
27 28 29 
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Again, from Eqs. (6.6) we have, 


F * 
^di 


T’ E> 

. + -qoi. di 

(x . -X* . ) (x . -X' . ) 
^ qi qi' ^ qi qi^ 


(6.64) 


Multiplying both sides of Eq. (6.64) by dEJj^/dt and integrating, 
we get the following equation after summing over all values of i 


m t 

- s / 


i=l t. 


■qi dt 


m 

y 

/ 

^0 


dF* 

°^di 

i=l 


dt 

m 

Y 

t 

/ 

^0 

t;.oi 


L* 

i=l 


Mt 


dt 


~ ^'^212"''^213 ’ view of Bqs.(6.42)-(6.43) 

(6.65) 

Thus, substituting from Eqs. (6.61), (6.63) and (6.65), the 
energy function in Eq. (6.29) can be written as 


W = Wj+W2i+W22+W23+Wi4+V<^5+Wi6 

where , 


I m n 

'^24 ^ I ^ ^^Gi-'^io^ ” -L ^%.i"%io^^ 

1=1 J -1 

m t dE* • 

’"25 = .=, { idi dt^ 

1=1 t^ 


'^26 = 


m 

- Z 
i=l 


t 

s 

t. 


A dt 

^qi dt 


( 6 . 66 ) 

(6.67) 

( 6 . 68 ) 


(6.69) 
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The expression for the energy function in Eg. (6.66) is much 
simple!' 3'^*^ requires less computations than that given in Eq.(6.29). 
W'^ accounts for the effect of flux decay and AVR and accounts 
for the effect of damper winding on the quadrature axis. It is 
obvious that both these terms would vanish when the generator is 
represented by the classical model and the expression, for the energy 
function will be the same as that used in Chapter 2. 

The SPEF defined in Eq. (6.66) can be evaluated using the 
computer program described in Chapter 2 with modifications to 
include the computation of the terras and in the sub- 

routine TEF for the computation of SPEF. 

6.4 NUMERICAL EXA^G^LES 
6.4.1 Description 

A 4-machine and a 10-machine system examples are considered 
here. The 10-machine system has already been described in previous 
chapters while the 4-machine system example is taken from [14]. The 
single-line diagram of 4-machine system is shown in Fig. K.l and 
the system parameters and the operating data are given in 
APPENDIX K. The disturbance considered is a three phase fault 
which is cleared followed by instantaneous reclosure of the line. 

The prefault and the post-fault configurations of the systems are, 
therefore, assumed to be the same. Impedance type load characteri- 
stics are assumed for all the busos. Critical energy is evaluated 
using the PEBS method. The following cases are investigated.: 
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Example 1 : 4-machine system 
Fault location : bus 3 

Aachines 2 and 3 which are close to the fault are modelled 
Ln detail considering voltage regulator, while machines 1 and 4 
are represented by classical models. 

Example 2 : 10-machine (NEW ENGLAND) test system 

Fault locations : i) bus 3 

i) Machine modelled in detail: 3 

ii) Machines modelled in detail: 1,3 

2 ) bus 4 

Machines modelled in detail : 4,5 

3 ) bus 8 

Machines 8 and 9 are represented with 
realistic models. 

The following is the data for AVR [71] : 

= 50 ; = 0.02 ; Kg = -0.037 

Tg = 0.146; Kg = 0.057; = 0.45 

0.015, = 0.6 ; 5.0 

Vin = 

For both the examples, the following cases are considered : 

A. Without Damper Winding 

1. without AVR ; a) with saliency; b) without saliency 

. a) with saliency; b) without saliency 


2. with AVR 
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B . With Damper Winding 

1, without AVR : a) with saliency, b) without saliency 
2» with AVR : a) with saliencyj b) without saliency 

6.4.2 Results and Discussion 
Critical Energies and Clearing Times 

Critical energy and critical clearing time obtained by 
prediction and digital simulation for various cases are 
presented in Table 6.1 for 4-machine system and Tables 6. 2-6 .4 
for 10-machine systems. The results in Tables 6. 3-6. 4 are 
obtained using the modified energy function neglecting the path 
dependent term in Vi22* energy component due to the active 
power of loads. 

From an examination of Tables, it is observed that the 
effect of flux decay alone is to reduce t^^ while the effect 
of transient saliency is to increase t^^. However, the effect 
of flux decay predominates the effect of transient saliency. 

The presence of the damper winding marginally improves the 
transient stability. The deleterious effect of flux decay 
is counteracted by the AVR. In the case of 4-machine system, 
tho critical clearing time is higher than that obtained using 
the classical model. This is not observed always in the case 
of 10-machine system. From Table 6.2 it is observed that 
providing AVR, on more than one machine is beneficial in 
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Table 6.3 Comparison of Critical Energy, W__ and Critical Clearing Time, t for lO-machine 
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increasing critical energies and critical clearing times. 

This is particularly true when transient saliency is considered. 
In general, it would be reasonably accurate to state that the 
classical model is a good and slightly conservative approxima- 
tion for the detailed model of a generator including AVR, 

In, practically, all the cases considered, the predicted 
values of critical clearing times are fairly close to the 
values obtained by digital simulation. The accuracy is slightly 
improved by including the kinetic energy correction as given in 
Chapter 3. The most significant observation is that the effects 
of flux decay, transient saliency, damper winding and AVR are 
correctly reflected in the predicted values of critical ener- 
gies and critical clearing times. This shows that the proposed 
SPEF can be used successfully in studying the stability chara- 
cteristics of a power system. 

Variations of Rotor Angles. Energy Components and other 
Variables 

Figs. 6. 2-6 .3 show the swing curves for example 1 
[B. 2— a)] for the stable and unstable cases for fault at bus 3 
of 4— machine system. Variations in the field voltage » 

Eq and E^ are given in Figs* 6. 4-6. 9 both for stable and un- 
stable cases. Figs. 6.10-6.19 show the variations in the total 
energy, kinetic energy, potential energy and its various compo- 
nents for the stable and unstable cases. 
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The variations in the rotor angles, El, energy 

and its various components for one case for fault at bus 3 
in the 10-machine system [B.2-a)] are shown in Figs. 6.20-6.37. 
Both stable and unstable cases are included here. 


In the case of the larger 10-machine system, when the 
system is unstable, there appear to be 3 distinct groups of 
generators that separate from each other (machine 2, machine 3 
and the remaining machines). In comparison with Fig. 3,11 
(Chapter 3), it appears that the provision of AVR on machine 1 
results in reducing the magnitude of oscillation within the 
coherent group. Even though it is not evident from the swing 
curves, the rotor oscillations tend to increase in magnitude 
for the critically cleared fault, which is reflected in the 
variation of or W2 (see Fig. 6.28). This problem could be 
due to dynamic instability and is normally cured by providing 
power system stabilisers which modulate the reference value 
(set point) of the AVRs. It is interesting to observe that 
this problem is not encountered in the case of 4-machine 
system. 

The action of AVRs is field forcing during the fault and 
immediately after the fault is cleared. This results in 
increase of E* . However, the AVR has little effect on EL which 

q 

tends to reduce from its initial value. But this has little 


effect on the stability. 
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Variations in and are shown alongwith the 

constituent compononts. has 5 components and their varia- 
tions appear to be correlated. has 2 components but their 

magnitudes are negligible compared to other components. This 
is not surprising as is due to the damper winding which 
has a very small effect on the transient stability. 


It is interesting to note that although the components of 
^25 large magnitudes (particularly, and W^g), the 

magnitude of is relatively small* This observation confirms 
the earlier statement that the use of the classical model is 
good approximation of the detailed model of the generator with 
AVR. This conclusion follows from the fact that the component 
^^25 vanishes when the classical model is considered. 


6.5 CONCLUSIONS 

In this chapter, a new SPEF is derived which takes into 
account the detailed models of generators including flux decay, 

it 

transient saliency, AVR and damper winding. This energy fun- 
ction is constant along' the post-fault trajectory. The poten- 
tial energy has 13 components which are finally shown equivalent 
to 6 components which can be speedily evaluated. A 4-machine 
and a 10-machine system examples are used to illustrate the 
application of the proposed SPEF for direct stability evaluation. 
The following are the conclusions of the analysis . 
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^cr with transient saliency, AVR and the 

inclusion of damper winding. 

ii) The predicted values of t^^ using the proposed energy 
function are close to the actual ones and the effects of 
flux decay » transient saliency» AVR and the damper winding 
can be accurately determined. 

iii) The variations of and have predominant effect on the 
system energy during the transient. 

iv) Variation of E^ and its effect on the transient energy is 
negligible. 

v) Although the detailed models of the generators including 
AVR are necessary for accurate evaluation of stability 
regions, the use of classical models is adequate as a first 
approximation. 


# 



CHAPTER 7 


CONCLUSIONS 


7.1 GENERAL 

Direct methods for transient stability analysis are 
potentially useful both as off-line tools for planning purposes 
as well as on-line dynamic security assessment. Although much 
research work has been reported in this area, the application 
of this method has been limited by conservativeness or accuracy 
of results, computational burden in determining stability re- 
gions and the use of simple (classical) models of generators and 
loads. Recent advances in the techniques of direct methods such 
as the use of energy functions in COI variables, the concept of 
controlling uep and PEBS in evaluating the region of stability 
have been directed at eliminating the limitations due to con- 
servativeness of results thereby enhancing the accuracy and 
reliability of these methods* However, the use of reduced net- 
works and classical machine models are major handicaps. 

Subsequent research work in the area has been aimed at 
removing the limitations due to system models. Structure 
preserving model which was first proposed by Bergen and Hill 
[ 52 ] is promising in that it allows detailed load representa- 
tions. Recently, Padiyar and Sastry [58] have developed a 
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'topological* enargy function which includes detailed generator 
and load representations . Preservation of network is also 
advantageous in on-line dynamic security assessment. However, 
the application of direct methods using structure preserving 
energy functions to on-r-line dynamic security assessment requires 
further developments due to i) more computational burden with 
nonlinear load models than with reduced models, ii) need to 
improve system models used and iii) the necessity to provide 
effective tools to the user for stability evaluation (both at 
the planning stage and in the operational level) of large 
systems. This thesis has been aimed at eliminating some of 
these limitations to make the direct methods acceptable to 
utilities. In this chapter, the major contributions of the 
thesis are reviewed and suggestions for further work are given. 

7.2 COMPUTATION OF SPEF IN DY^IA^UC SIMULATION PROGRAM 

Although efficient computer programs to perform digital 
simulation for stability evaluation are available and capable 
of accommodating detailed generator models, the non availabl 1 y 
of effective programs for the computation of energy functions 
has prevented extensive application of direct methods to system 
planning. However, the existing programs. If suitably modified, 
can be used to study the transient behaviour of power systems 
through direct methods or simulation. In Chapter 2, an exist- 
ing digital simulation program [66] has been i) modified to 
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include voltage dependant loads and then ii) augmented to 
include the computation of SPEF and its- various components by 
incorporating two new subroutines PFSOL and TEF. While sub- 
routine- PFSOL solves the post-fault network equations at the 
end of each stop of integration of the swing equations des- 
cribing the system dynamics, subroutine TEF evaluates the SPEF 
and its components utilising the bus voltage magnitudes and 
the bus angles obtained as outputs of PFSOL. The resultant 
program can, thus, be conveniently used to analyse the tran- 
sient stability of power systems both by digital simulation and 
by direct methods based on the computation of SPEF. 

The analysis of the variations of kinetic and potential 
energy components is helpful in deciding system stability 
thereby avoiding the need to scrutinise the voluminous results 
in the form of swiing curves. 

7.3 STRUCTURE PRESERVING ENERGY FUNCTIONS WITH KNOWN MODES 

OF INSTABILITY 

In the application of PEBS method, it is assumed that 
the kinetic energy gained by the faulted system at the instant 
of clearing gets fully converted into the potential, energy of 
the post— fault system. Though the assumption is true for a 
two-machine system, it has been shown [59] that the total 
kinetic energy need not be converted into potential energy in 
the case of multi-machine systems. It is the kinetic energy 
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at clearing corresponding to a particular step-out mode that is 
responsible for system separation [32], in this thesis, a novel 
structuiT© piTGSGjrvinQf function hss boon dovslopod (in 

Chapter 3) using CX)I variables under the assumption that the 
system splits into r groups (r 2 2) on the onset of instability 
and that the generators within a group swing together. The 
kinetic energy component of the energy function excludes the 
energy due to inter-generator swings within the group and thus 
gives better estimate of the energy associated with system 
separation. Although this energy component can be derived from 
the SPEF used in Chapter 2 [54] , the proposed SPEF is derived 
directly from a system model based on the mode of instability. 

Numerical examples demonstrate that the kinetic energy 
correction by the application of the proposed SPEF inproves 
the predicted results only marginally. This is because the 
kinetic energy component relative motion of 

rotors of generators in a coherent group is not zero and does 
not remain constant along the post-fault trajectory. There is 
an exchange of energy between the two components of the kinetic 
energy term. This fact seems to be ignored in the previous 
literature on the subject. 

The presence of non-constant active power loads results in 
path dependant terms in the energy function. These terms can be 
integrated numerically without difficulties. However, it is 
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angles become implicit functions of the faulted generator rotor 
angle* While the trajectory of the faulted generator rotor 
which accelerates monotonically during the fault-on period, can 
be determined by increasing its angle from its initial value by 
an arbitrarily selected step size, the remaining system varia- 
bles required for the computation of SPEF can be obtained by 
solving a set of algebraic equations* Further, the assumption 
of constant accelerating power for the faulted generator till 
the instant of fault clearing, speeds up the computation of. 
energy at clearing thereby reducing the overall computation 
time for TEM. The algorithm for the calculation of TEM is 
validated by the correct prediction of stability or instability 
of the systems investigated and the ranking of the contingencies 
considered* 

The work reported in this thesis is first step in the work 
required to implement the application of SPEF for on-line dyna- 
mic security assessment* The elimination of the computation of 
the faulted trajectory in detail is significant not only in 
terms of reducing the computations but allowing the introduction 
of sensitivity analysis for handling changes in loading condi- 
tions and the network configurations. 

The algorithm given in Chapter 4 is based on the assumption 
that the faulted generator tends to separate from the rest. As 
seen from the results in Chapters 2 and 3, this assumption is 
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valid for most of the cases. Where it is not applicable (for 
example, fault at the terminals of generator 10 of the 10— 
machine system), the contingency would be, usually, not criti- 
cal and the results of the security assessment would not be 
affected. 

7.5 STRUCTURE PRESERVING ENERGY FUNCTIONS WITH DETAILED 

SYSTEM MODELS 

A novel method has been presented in Chapter 5 to account 
for transmission line resistances in the SPEF. The method 
exploits the equivalence of a lossy network having same G/b 
ratio for all its elements to a lossless network with a new 
set of power injections. The new SPEF developed using 001 
variables has a path dependent term which can be numerically 
integrated without affecting the accuracy in the predicted 
results. The error due to the approximation in assuming con- 
stant G/B ratio for all the lines is found to be small. 

Reduced networks where the elimination of load buses 
introduces transfer conductances which are predominant, require 
the inclusion of path dependent term in the energy function. As 
remarked earlier, the presence of largo path dependent terms is 
not desirable and perhaps explains the reason for the failure 
of the PEBS method in accurately predicting the critical energy 
sometimes as reported in the literature. The accuracy of the 
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SPEF developed in Chapter 5 is very good perhaps due to 
relatively small path dependent term in the energy function. 

The effect of line resistances seem to improve the critical 
clearing time and this is accurately predicted by the SPEF. 

A SPEF has been developed in this thesis considering flux 
decay » transient saliency^ exciter j AVR and a damper winding 
on the quadrature axis. Recently, a 'topological* energy fun- 
ction (TEF) considering realistic generator models (these do 
not include the damper winding) has been reported in ref, [58], 
However, the TEF in ref. [58] has a time derivative which is 
negative semi-definite along the post-fault trajectory and 
gives conservative estimates of the region of stability. The 
SPEF proposed in this work removes the conservativeness of the 
results* Moreover, a general expression for the SPEF is deri- 
ved which simplifies the computation of the energy function. 

From a study of two system examples considered, the following 
conclusions can be derived ; 

i) The effect of transient saliency, AVR and damper winding is 
to increase the stability region. However, the effect of 
damper winding on system stability is not very perceptible. 

ii) The use of SPEF gives reliable results and correctly shows 
the effects with control equipments, 

iii) For stability analysis, the use of a classical model for 
generators is a good approximation for the detailed models. 
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7.6 SUGGESTIONS FOR FURTHER VJORK 

Based on the observations made in the previous sections, 
the following suggestions are made for further work in this 
area : 

i) An algorithm for accurately predicting the mode of insta- 
bility for a given contingency without computing the system 
trajectory is not yet available and requires investigation. 

ii) In this work, a beginning has been made in the application 
of SPEF for on-line dynamic security assessment. This work 
can be extended to further simplify the computations and 
implement the use of SPEF in energy control centres. This 
would require the following developments : 

a) The use of sensitivity analysis for changing load 
conditions and the network configurations. 

b) At present, the direct method considers only a single 
element such as fault followed by its clearing. An 
emergency state can lead to a series of events. It is 
worthwhile examining whether SPEF can handle the 
analysis of such events. 

c) Before actual implementation extensive testing on large 
systems is required. 
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APPENDIX A 

DERIVATION OF SYNCHRONOUS GENERATOR EQUATIONS [2, 66] 
A.l PARK EQUATIONS 


The Park equations of synchronous generator with a d-axis 
field winding and a q-axis damper winding are given below : 

Direct axis : 


*cl - ^ad^f ^d^d 

(A.l) 


(A.2) 

d* . 

''d = ■ dt" - Vd 

(A.3) 

dUf 

''f = ar + 

(A.4) 

Quadrature axis : 


^q “ “^aq^kq ^q^q 

(A. 5) 

^kq ~ ^kkq^kq ” ^aq^q 

(A. 6) 

d’J?^ 

q dt a q o 

(A.7) 

0 = cTT^ ^kq^kq 

(A.8) 
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A. 2 STATOR EQUATIONS 

In this section the stator equations which are suited 
for use in the generator model are derived. For the purpose 
of transient and dynamic stability studies the stator trans- 
former voltages can be neglected. Eqs. (A. 3) and (A. 7), 
therefore, become : 




- a)W 


(A.9) 


V = -r i + (0^; , 
q a q u 


(A.IO) 


The q-axis stator equation is firstly derived by eliminating 
i_ from (A.l) and (A. 2), 

I 

Vad^^ff ‘^^^d'^ad'^^ff^^d (A.H) 

Define the d-axis transient inductance. 


- L^d/L.f . 

and the q-axis transient voltage, 


(A.12) 


and substitute into (a*11)» 

+ mLVji^ 

Substitute (A. 14) into (A.IO), 

Vq = Eq "■ ^a^q 


(A.13) 


(A.i4) 


(A.15) 
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The d—axis stator equation is now derived by firstly 
ing ij^q from (A. 5) and (a«6), 


“*c, “ *kq‘-aq/4kq + ““•q-‘'V''ttq)iq 
Define the q-axis transient inductance, 

^q " ^q “ ^aq/^kkq 
and the d-axis transient voltage, 


^d ^kq ^aq'^^kkq 

and substitute into (A. 16), 

wf = -E’, + wL'i^ 
q d q q 

Substitute (A.19) into (A.9), 

^d ~ ^d ^^q^q ” 

A. 3 ROTOR WINDING EQUATIONS 


The d-axis field winding equation is derived as 
Multiply (A.4) by coL^^/r^ 
coL. 


V. 


■ad . “^ad '''*f 


Define , 


0)1 


h = 


ad 


‘^^ad 
Ic 2? f 

Kfr^ 


eliminat- 

( A.16) 

(A.17) 

(A.18) 

(A.19) 

(A.20) 

follows: 

(A.21) 


V 


(A.22) 
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From ( A*2) , 


if ~ ^ad^d^^ff 


(A.23) 


Substitute (A. 22) and (a. 23) into (A. 21), 


kfEf 


(oL , wL' 

ad m 

f L 


ff d Tf dt 


(A.24) 


Define the d-axis open circuit transient time constant, 
Tdo = 


ff' -"f 


Combine (A. 12), (a» 13), (A.25) into (A.24), 

dE' 

kfEf = E' - coCLj-L^jlf, + TJi, 3^ 


(A.25) 


(A.26) 


The q-axis damper equation is derived as follows : Multiply 
(A.8) by 


0 = .L 1 . 

‘-aq kq dt 

From (8.6) , 

^kq “ '^kq'^^kkq ^aq^q^^kkq 


(A.27) 


(A.28) 


Substitute (A.28) into (A.27), 

2 wL d f 


wL ojL 

0 = r— ^ . + r— 

^kkq ^kkq 


ac[ . , _aa — isa 

^q ^ r,. . dt 


^kq 


(A.29) 


Define the q—axis open circuit transient tim© constant. 


^qo ^kkq'^^kq 


(A.30) 


Combine (a.18) , (A. 19) and (A.30) into (A.29), 


dE' 


0 = E' + a)(L -L* )i„ + T' 


q q 9 


■qo dt 


(A.31) 
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A. 4 AIR GAP VOLTAGE EQUATIONS 

The d and q-axis air gap voltages are defined by, 

^aq = '"*ad = “(^d + (A.32) 

and + 41 ^) U.33) 

The equations used in the generator model are obtained as 
follows : Substitute (A. 14) into (a. 32), 

E = E’ + w(L' + L.)i, (A.34) 

aq q did 

Substitute (Aal9) into (A. 33), 

^ad = 

In Sections A»1 to A»4 the generator model equations have been 
developed in terms of inductances* For transient stability 
studies the deviations in frequency are small and one can sub- 
stitute a)L = X ,(x is the reactance at synchronous frequency) 
in the equations derived. 

A. 5 STEADY STATE GENERATOR EQUATIONS 

The vector diagram for synchronous generator in the steady 
state is shown in Fig. A.l. The voltage E^ is defined as 




(A. 36) 
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During the steady state lies on the quadrature axis, since 


Im[E 1 = + r^I^ + X I 

H q q 

(A.37) 

and from (A. 5) and (A. 9) 


d ad q q aq kq 

(A.38) 


In the steady-state, =0 and, therefore, Im[E ] = 0. This 
property of E is used in calculating the generator initial 

Hi 

conditions since the terminal voltage and current are given 
in terms of the synchronous reference frame. The position of 
the generator q-axis is given by calculating, 



(r +jx )T 


(A.39) 


and 6 = Arg(E^) 


(A.40) 


From the vector diagram it can also be seen that in the steady 
state , 

= ^q - (’‘q-='d>^d 

The other equations which are used in calculating the generator 
initial conditions are obtained by setting derivatives to zero 


in (A»26) and (A*31) ^ 

Ef = [Eq “ (^cJ”^d^^d^'^^f 


^d 




(A.42) 


(A. 43) 
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APPENDIX B 


network solution with nonlinear voltage dependent loads [69] 


In the presence of nonlinear voltage dependent loads, the 
computation of faulted trajectory requires iterative solution 
of the nonlinear algehraro network equations* An algorithm 
giving fast convergence is discussed here. 


The admittance matrix [Y^^] for the system network is 
formulated (excluding machine reactance). A generator is 
represented by its Norton's equivalent i.e., a current source 
in parallel with a fictitious admittance to be obtained as 


fict 

^i 


”2 


ai 


+ x 


di"qi 


(B.l) 


A load is represented as a current source in parallel with a 
constant admittance (see Fig. B.l). Admittance, y^j at bus j 
is evaluated from the initial load flow and hence, initially, 
I^j is zero. The admittance matrix is augmented to 

include at the generator bus i and y^j at the load bus j . 

Let the augmented admittance matrix be [Y]. Then the network 
equations are : 



ill) = Cy] [y] 


where I 


Gi 


_ c^ict 
- 


fict 
^i ’ 


= 1 , 2 . 


• • • J 


(B.2) 


(B«3) 
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'Fi. c 

e!" is the voltage behind the fictitious admittance yf of 
ith generator and is calculated from 


1 


=fict 


where 


r, ^di^ /■=.* — -ioo 

1 " ; ' "3 1,. . W.2 


"i = ^ii 


+ , i = 1,2, ..., m 


m 

(B.4) 

(B.5) 


vector [lj^3 represents the current injections due to loads. 
Current injection 1^^^ at a load bus j can be obtained from 



(B.6) 


The admittance matrix changes whenever the network changes 
due to fault, line switching, etc., otherwise it remains con- 
stant at all iteration steps. 

If all the loads are represented as impedances, the Eg. 
(B.2) would be non-iterative with [l^^] =0. If the loads are 
not modelled as impedances, then [l.^] is adjusted iteratively 
to account for the difference between the quadratic voltage 
dependence of the impedance and the specified voltage dependenct 
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APPENDIX C 


TIME DERIVATIVE OF ENERGY FUNCTION WITH NONCONSTANT MECHANICAL 

POWER 


The structure preserving energy function including the 
nonconstant mechanical power can be written as 


W = + W 2 = + S W2i 


(C.l) 


where 


l"" 2 

W, = w 2 M.o)^ 

i=l ^ ^ 

m t de. 

»21 = - 3r 

o 

. n t d9. 

^22 •-! { ^Pi^^i^ ^ 

i-x 

n "^i f 4(V.) 

^23 = ^ ^VT^ "^^i 

23 i=i V.^ ^i 


(C.2) 


(C.3) 


(C.4) 


(C.5) 


2 [E?+V^-2E^V^ cos(0^~<P^) - CE?+V?q 


-2E.v.„ cosCe.„-<P.,))] 5 ^ 


(C.6) 


n n 




<P^ -! J 


(C.7) 
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The time derivative of the function W in Eq. (C 4 I) can be ♦ 
shown to be zero by direct verification as follows : 

Differentiating with respect to we have, 
dW, 

^cirr ~ » i='l,2,«.»,m (C«£ 

Taking the partial derivative of W 2 with respect to ^^*^^** 
and t respectively, we get. 


dW^ 


E.V- sin(e.-<p. ) 

1 = 1,2, 


X 


di 


ei 


• ♦ > 


m 


(C.9) 


dW, 


3v. 


2 _ 


1 ^ 1 . - cos(ej_-<p.) 


^di 


^di 


J 

r ~ 1 , 2 , m 


n , 

= - 2 B..V.V. cos (P. ^ , i =m+l,m+2,... 

-L i=l -^J J '' 3 ^ 


~ vT" ^^i ^ •••} 


= 0 


(C.IO) 

(C.IO) 


dW^ E.V. sin(0.-(P. ) n 

L_i — , i — 1- + 2 B. .V.V. sin (P, . , i=l,2, 


d<P 


X * 

^di 


_i 13 1 J 


3=1 


13 


, m 


n 


= Z B. .V. V . sin (p . . , i = m+1 , m+2 , . . . , n 
13 1 3 13 


= Pi» 


1 ^ f'2. f • • • f n 


(c.ll) 
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dW 2 m d 0 . n d<J). 

^fr + ^Li Ht“ 


(C. 12 ) 


From Eqs. {C. 8 ), (C. 9 ) , (C.ll) and (C. 12 ), we get after 
rearranging the terras, 

m dw, cb. m aw^ de. n dViL d<p. dW^ 

SS 7 dT 5 ^ ar 357 Hr sr 

m . j. n d<p. m n d<P^ 

= _! Mj^Vi + \ + .S, Pi dT - ■" ,i, '’Ll Hr 

1—1 1=1 1=1 1—1 J- 


m M.u). n dg). 

= " 1=1 TiT ‘’®' 1=1 ^'’i'’'''Li’ Hr 


Pcoi “ 


M 


T i 


m n d<p. 

E E (p +P ) ^ 

1 -i 1=1 


= 0 
Hence, 


if = Ht ("rt) 


m aw, do). .dWn de. n dW^ d<P. dW2 

C (jr af - 3^ Hr) h 57 HT + sr^ 


n aw. dv. 


(C. 13 ) 


= 0 , in view of Eqs. (C.IO) and (C. 13 ) 


(C. 14 ) 
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APPENDIX D 

PEBS METHOD OF DETEPJ^AINING STABILITY REGION 

As discussed earlier in Chapter 1, the potential energy 
boundary surface is the surface in the angle space of the points 
corresponding to the first maxima of the potential energy with 
respect to the stable equilibrium point where its directional 
derivative is zero. The method of computing stability region 
using PEBS is explained below : 

The basic idea of PEBS method is to replace the study of 
the following model of power system consisting of m generators 

9 f i = 1^2} ni (D.la) 

, M. 

^i'^i ~ ^mi”^ei"' ^COI* i - 1,2, m (D.lb) 

by the study of simple gradient system [77] 

dWo 

S {D.2) 

d © 

In the PEBS method we calculate the critical energy from the 
system (D.2) rather than from system (D.l) using the properties 
of the gradient system. However, it has been shown in ref .[77] 
that the boundary of the stability region of the system descri- 
bed by Eqs. (D.l) is related to that of the gradient system of 
Eq. (D.2). 
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The structure preserving energy function can be written as 


W(g.i 

>1 

31 

, <Lf 

t) = Wj_(w) + W, 

2(§.. 

V, 

t) 



(D.3) 

Hence, 










dWi 

dt 

m 

= 2 
i=l 


do)^ 

dt 







(D.4) 

ciW2 

dt 

m 

= 2 
i=l 

dW2 

69-^ 

d0. 
dT "• 

n dW« 

1=1 H 

d^i 

dt 

aw2 

dV. 

dW, 

3t" 




m 

= 2 
i=l 

dV«2 

60^ 

de. 

1 . 

n awo 

1=1 1 

d<p^ 

dt 

' aw^ 

dt J 

m 

= 2 
i=l 

SWg 

357 

de. 











(D.5) 

dW 

dt 

m 

= E 

i=l 

M^oo. 

+ . 

a 

m aw^ . 

^ 3r ®: 
.=1 1 

L 

m 

2 (M^ 

i=l 


dW, 

357 ' “i 

1 

(D.6) 

Let 

vector S 

== 

L ^2^2 

• • • 

M S) ]■ 
m m-* 

r 





dW2 

dWo 

aw^ 


awo T 





and 


= [557 


• 

m 






Therefore , 


dt 


aw2 


= 0 


(D.7) 


for a conservative system. From this it follows then, 

= - W 


(D.8) 
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The direction of the torque vector S is in a direction 
negative to the gradient of W 2 » i.e., S = - 

By definition VW 2 is in a direction normal to W 2 = 
constant in the angle space. Hence, torque is always ortho- 
gonal to the equipotential surfaces and in the direction of 
decreasing W 2 * Consider a trajectory crossing the PEBS. 

Inside the region surrounded by PEBS and containing the sep, 
the projection of S on the velocity vector w is in a direction 
opposing the motion, i.e., towards sep. Outside PEBS, proje- 
ction of S is in the direction of the velocity vector u. Thus, 
the dot product of VW 2 with velocity vector to changes sign in 
the vicinity of the PEBS. 

From Eq. (D.7) , we get, 

dVL m d0. 

s-a = - • a = - .z 35^ cTT = 

1 =1 1 

in view of Eq. (D.5) . 

Hence, W 2 > 0 inside the PEBS and < 0 outside it and vice versa 

o 

is true with Vi/^. This implies that inside the stable region the 
potential energy increases and the kinetic energy decreases. The 

o « 

opposite is true outside the PEBS. Consequently, when W 2 or 
changes sign can be interpreted to be the instant when the 
trajectory crosses the PEBS [29]. 
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The computation of critical energy, W__ and critical 

C X 

clearing time, t^^ by the PEBS method consists of the follovd.ng 
procedure : 

1) Compute the fault-on trajectory. 

2) Monitor , the maximum value of on the fault-on 

trajectory. Critical energy, V/ = . 

OX 

3) Find t^^^ such that 

W(t^l) = V«i(t^i) + W2(tel) = Wor 
The t^j^ of Step 3 is the PEBS method solution for the 
critical clearing time. 


4) 
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APPENDIX E 


DATA FOR THE 7-MAGHINE (CIGRE) TEST SYSTHA 


Number of machines 
Number of buses 
Number of lines 
Number of loads 
Base MVA 


Base kV : 225 


Table E.l Machine Constants 



7 


0 .0568 


0.087 




282 


Table E.2 Line Data 


Line 

No. 

From Bus 

To Bus 

Line reactance 
in ohms 

Half-line 

charging 

(i^s) 

1 

2 

3 

24.5 

200 

2 

2 

4 

24.5 

;oo 

3 

3 

5 

62.6 

200 

4 

5 

10 

32.3 

300 

5 

3 

4 

39.5 

300 

6 

3 

9 

28.0 

200 

7 

4 

1 

10.0 

200 

8 

4 

6 

24.8 

100 

9 

4 

9 

97.0 

200 

10 

4 

10 

33.0 

300 

11 

6 

8 

31.8 

200 

12 

7 

8 

39.5 

300 

13 

8 

9 

97.0 

200 
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Table Ei3 Operating Data 


Bus 

No, 

Voltage 
Magnitude 
in pu 

Bus angle 
in deg 

Generation 

Load 


Active 
pov^^er 
in pu 

Reactive 
power 
in pu 

Ac tive 
power 
in pu 

Reactive 
power 
in pu 

1 

1 .0343 

-1 .4640 

' 2.30 

1.0 

0.0 

0.00 

2 

1 .0448 

-0 .1852 

2.00 

0.7346 

0.0 

0.00 

3 

1 .0363 

-1.6300 

2.56 

0.8536 

0.0 

0.00 

4 

1 .0141 

-3.9461 

3.00 

1.7880 

6.5 

4.05 

5 

1 .0257 

-7 .0880 

1.20 

1 .0602 

2.0 

1.20 

5 

1 .0125 

-2.8794 

1.60 

0.2919 

0.8 

0.30 

7 

1 .0103 

-0.6210 

1.74 

0.3832 

0.9 

0.40 

8 

1 .0017 

-4.3339 

0.00 

0.0 

1.0 

0.50 

9 

0.9829 

-7 .12 

0.00 

0.0 

2.3 

1.40 

10 

1 .0148 

-7.15 

0.00 

0.0 

0.9 

0.45 
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APPENDIX F 

DATA FOR THE 10-MACHINE (NEW ENGLAND) TEST SYSTEM 


Number of machines : 10 
Number of buses : 39 
Number of lines ; 34 
Number of transformers: 12 
Number of loads : 19 


Table F.l Machine Constants 


Generator 

No, 

Bus 

No. 

H 

in sec 

• in pu 

in py 

^d 

in pu 

in sec 

1 

1 

30.3 

0,0647 

0.1543 

0.1812 

7.5 

2 

2 

500.0 

0.006 

- 

- 

«■» 

3 

3 

35.8 

0.0531 

0.1261 

0.1493 

5.8 

4 

4 

23.6 

0 .0436 

0.1040 

0.1224 

5.1 

5 

5 

26.0 

0.1320 

0.3142 

0.3696 

10.0 

6 

6 

34.8 

0.05 

- 

~ 


7 

7 

26,4 

0.049 

- 

- 


8 

8 

24.3 

0.057 

0.1359 

0.1596 

6.0 

9 

9 

34.5 

0.057 

0.1359 

0.1596 

6.0 

10 

10 

42.0 

0.031 

•• 
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Table F.2 Line Data 


Line 

No. 

From 

Bus 

To 

Bus 

Line 

impedance 

Half-line 
Charging 
in pu 

Resistance Reactance 
in pu in pu 

1 

2 

3 

4 

5 

6 

1 

30 

31 

0.0035 

0 .0411 

0.3494 

2 

30 

2 

0 .0010 

0 .0250 

0.3750 

3 

31 

32 

0.0013 

0 .0151 

0.1286 

4 

31 

25 

0 .0070 

0 .0086 

0 .0730 

5 

32 

35 

0.0013 

0.0213 

0.1107 

6 

32 

18 

0 .0011 

0.0133 

0.1069 

7 

33 

34 

0 .0008 

0.0128 

0.0671 

8 

33 

14 

0 .0008 

0 .0129 

0.0691 

9 

34 

35 

0 .0002 

0 .0026 

0.0217 

10 

34 

37 

0 .0008 

0.0112 

0 .0738 

11 

35 

36 

0 .0006 

0 .0092 

0.0565 

12 

35 

11 

0 .0007 

0 .0082 

0.0695 

13 

36 

37 

0 .0004 

0 .0046 

0.0390 

14 

37 

38 

0 .0023 

0.0363 

0.1902 

15 

38 

2 

0 .0010 

0 .0250 

0.6000 

16 

39 

11 

0 .0004 

0 .0043 

0.0365 

17 

39 

13 

0 .0004 

0.0043 

0.0365 

18 

13 

14 

0 .0009 

0.0101 

0.0862 

19 

14 

15 

0 .0018 

Q.0217 

0.1830 

20 

15 

16 

0 .0009 

0 .0094 

0.0855 


contd 
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t 5 


1 

2 

3 

4 

5 

6 

21 

16 

17 

0.0007 

0 .0089 

0 .0671 

22 

16 

19 

0.0016 

0 .0195 

0.1520 

23 

16 

21 

0.0008 

0.0135 

0.1274 

24 

16 

24 

0.0003 

0 .0059 

0.0340 . 

25 

17 

18 

0.0007 

0.0082 

0 .0659 

26 

17 . 

27 

0.0013 

0.0173 

0.1608 

27 

21 

22 

0.0008 

0.0135 

0 .1274 

28 

22 

23 

0.0006 

0 .0096 

0.0923 ■ 

29 

23 

24 

0.0022 

0.0350 

0.1805 

30 

25 

26 

0.0032 

0 .0323 

0.2565 

31 

26 

27 

0 .0014 

0 .0147 

0.1198 

32 

26 

23 

0 .0043 

0 .0474 

0.3901 

33 

26 

29 

0.0057 

0.0625 

0.5145 

34 

28 

29 

0 .0014 

0 .0151 

0.1245 
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Table F.3 Transformer Data 


Transformer 

No. 

From 

Bus 

To 

Bus 

Impedance 

Turns 

ratio 

Resistance 
in pu 

Reactance 
in pu 

1 

12 

11 

0 -0016 

0.0435 

1.006 

2 

12 

13 

0.0016 

0.0435 

1.006 

3 

35 

1 

0.0 

0 .0250 

1.070 

4 

39 

3 

0.0 

0 .0200 

1 .070 

5 

19 

4 

0 .0007 

0 .0142 

1 .070 

6 

20 

5 

0 .0009 

0.0180 

1.009 

7 

22 

6 

0.0 

0.0143 

1.025 

8 

23 

7 

0 .0005 

0.0272 

1.0 

9 

25 

8 

0 .0006 

0.0232 

1.025 

10 

31 

10 

0.0 

0.0181 

1.025 

11 

29 

9 

0 .0008 

0,0156 

1.025 

12 

19 

20 

0 .0007 

0.0138 

1.060 
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Table F.4 Operating Data (Including line resistance) 


Bus 

No. 

Voltage 
Magnitude 
in pu 

Bus 
angle 
in deg 

Generation 


Loads 

Active 
power 
in pu 

Reactive 
power 
in pu 

Active 
power 
in pu 

Reactive 
power 
in pu 

1 

2 

3 

4 

5 

6 

7 

1 

0 .9820 

0.0 

5.484 

2 .157 

0 .092 

0.046 

2 

1 .0300 

-10.6620 

10.000 

0.983 

11 .040 

2.500 

3 

0 .9831 

2.0860 

6 .500 

2 .230 

0.0 

0.0 

4 

0 .9972 

3 .4710 

6 .320 

1.444 

0.0 

0.0 

5 

1 .0123 

2.4320 

5 .080 

1.816 

0.0 

0.0 

6 

1 .0493 

4.6790 

6.500 

2.596 

0.0 

0.0 

7 

1 .0635 

7 .3490 

5.600 

1.323 

0.0 

0.0 

8 

1 .0278 

1 .7331 

5 .400 

0.1184 

0.0 

0.0 

9 

1 .0265 

7 .1022 

3.300 

0.3486 

0.0 

0.0 

10 

1 .0475 

-4 .0150 

2.500 

1.6155 

0.0 

0,0 

11 

1 .0089 

-6 .7821 

0.0 

0.0 

0.0 

0.0 

12 

0.9958 

-6.7621 

0.0 

0.0 

0.075 

0.880 

13 

1 .0097 

-6.6363 

0.0 

0.0 

0.0 

0.0 

14 

1 .0053 

-8 .2621 

0.0 

0.0 

0.0 

0.0 

15 

1 .0033 

-8 .5024 

0.0 

0.0 

3.200 

1.530 

16 

1 .0175 

-6 .9806 

0.0 

0.0 

3.294 

0.323 

17 

1 .0228 

-8 .0752 

0.0 

0.0 

0.0 

0.0 

18 

1 .0220 

-8 .9790 

0.0 

0.0 

1.580 

0.300 

19 

1 .0446 

-1 .7585 

0.0 

0.0 

0.0 

0.0 

20 

0.9883 

-2.7650 

0.0 

0.0 

6.280 

1.030 


contd 
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21 

1 .0215 

- 4.5945 

0.0 

0.0 

2.740 

1.150 

22 

1 .0432 

- 0.3143 

0.0 

0.0 

0.0 

0.0 

23 

1 .0369 

- 0 .5562 

0.0 

0.0 

2.745 

0.846 

24 

1 .0149 

- 6 .8910 

0.0 

0.0 

3.086 

0.922 

25 

1 .0549 

- 5 .0647 

0.0 

0.0 

2.240 

0.472 

26 

1 .0468 

- 6 .2636 

0.0 

0.0 

1 .390 

0.170 

27 

1 .0298 

- 8.2627 

0.0 

0.0 

2.810 

0.755 

28 

1 .0474 

- 2.7362 

0.0 

0.0 

2.060 

0.276 

29 

1 .0481 

0.0313 

0.0 

0.0 

2.835 

0.269 

30 

1 .0464 

- 9 .0820 

0.0 

0.0 

0.0 

0.0 

31 

1 .0460 

- 6.4410 

0.0 

0.0 

0.0 

0.0 

32 

1 .0243 

- 9.2300 

0.0 

0.0 

3.220 

0.024 

33 

0.9985 

- 10 .2040 

0.0 

0.6 

5.00 

1.840 

34 

1 .0014 

- 9 .0998 

0.0 

0.0 

0.0 

0.0 

35 

1 .0040 

- 8 .4120 

0.0 

0.0 

0.0 

0.0 

36 

0.9934 

- 10 .6170 

0.0 

0.0 

2.338 

0.840 

37 

0 .9925 

- 11.1200 

0.0 

0.0 

5.220 

1.760 

38 

1 .0268 

- 10 .8900 

0.0 

0.0 

0.0 

0.0 

39 

1 .0133 

- 5.9400 

0.0 

0.0 

0.0 

0.0 
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Table F.5 Operating Data (neglecting line resistance) 


Bus 

NO . 

Voltage 
Magnitude 
in pu 

Bus angle 
in deg 

Generation 


Loads 

Active 
power 
in pu 

Reactive 
power 
in pu 

Active 
power 
in pu 

Reactive 
power 
in pu 

i 

2 

3 

4 

5 

6 

7 

1 

0.982 

0 *0 y 

5.045 

1.9164 

0 .092 

0.046 ■ 

2 

1.030 

-9.538 

10.0 

0.684 

11.04 

2.5 

3 

0.9831 

2.9 

6.5 

2.159 

0.0 

0.0 

4 

0.9972 

4.806 

6.32 

1.57 

0.0 

0.0 

5 

1 .0123 

3 .834 

5.08 

1.852 

0.0 

0.0 

6 

1.0493 

5.903 

6.5 

2.547 

0.0 

0.0 

7 

1.0635 

3.584 

5.6 

1 .376 

0.0 

0.0 

8 

1 .0278 

2 .853 

5.4 

0.46 

0.0 

0.0 

9 

1 .0265 

8.378 

8.3 

0.597' 

0.0 

0.0 

10 

1 .0475 

-2.7 

2.5 

1.256 

0.0 

0.0 

11 

1.0116 

-5.99 

0.0 

0.0 

0.0 

0.0 

12 

0.9989 

-5.98 

0.0 

0.0 

0.075 

0.88 

13 

1 .0129 

-5.79 

0.0 

0.0 

0.0 

0.0 

14 

1 .0122 

-7.34 

0.0 

0.0 

0.0 

0.0 

15 

1.0131 

-7.42 

0.0 

0.0 

3.2 

1.53 

16 

1 .0256 

-5.8 

0.0 

0.0 

3.294 

0.323 

17 

1 ,0320 

-6.87 

0.0 

0.0 

0.0 

0.0 

10 

1 .0322 

-7.774 

0.0 

0.0 

1.58 

0.30 

19 

1 .0475 

-0.468 

0.0 

0.0 

0.0 

0.0 

20 

0.9924 

-1.436 

0.0 

0.0 

6.28 

1.03 
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1 

2 

3 

4 

5 

6 

7 

21 

1 .0270 

- 3.443 

0.0 

. 0.0 

2.74 

1.15 

22 

1 .0439 

0.913 

0.0 

0.0 

0.0 

0.0 

23 

1 .0382 

0.655 

0.0 

0.0 

2.745 

0.846 

24 

1 .0230 

- 5 .71 

0.0 

0.0 

3 .086 

0.922 

25 

1 .0503 

- 3 .978 

0.0 

0.0 

2.24 

0.472 

26 

1 .0526 

- 5 .01 

0.0 

0.0 

1.39 

0.17 

27 

1 .0393 

- 7 .037 

0.0 

0.0 

2.81 

0.755 

28 

1 .0530 

- 1 .44 

c.o 

0.0 

2.06 

0.276 

29 

1 .0508 

1.31 

0.0 

0.0 

2.835 

0.269 

30 

1 .0496 

- 7.842 

0.0 

0.0 

0.0 

0.0 

31 

1.0524 

- 5.112 

0.0 

0.0 • 

0.0 

0.0 

32 

1 .0343 

- 8.1 

0.0 

0.0 

3.22 

0.024 

33 

1.0073 

- 9.286 

0.0 

0.0 

5.0 

1.84 

34 

1 .0074 

- 8.34 

0.0 

0.0 

0.0 

0.0 

35 

1 .0089 

- 7.686 

0.0 

0.0 

0.0 

0.0 

36 

1 ,0006 

- 9.863 

0.0 

0.0 

2.338 

0.84 

37 

1 .0003 

- 10.345 

0.0 

0.0 

5.22 

1.76 

33 

1 .0299 

- 9.861 

0.0 

0.0 

0.0 

0.0 

39 

1 .0148 

- 5.112 

0.0 

0.0 

0.0 

0.0 
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APPENDIX G 

DATA FOR THE 13-MACHINE SYSTEM 


Number of 

buses 

= 71 

Number 

of 

machines 

= 13 

Number of 

lines 

= 94 

Number 

of 

shunt capacitors 

= 23 

Base MVA 


= 200 






Table ,G.l : Machine constants. 


Gen. 

No. 

Bus 
No . 

Rated 

MVA 

H 

(in sec) 

T' 

Mo 

(in sec) 

(in pu) 

^d 

(in pu) 

1 

1 

1175.0 

12.9 

7.0 

0.36 

0.0488 

2 

3 

665.0 

9.0 

6.5 

0.535 

0.069 

3 

6 

110.0 

1.923 

6.0 

1.638 

0.5278 

4 

8 

335.0 

6.648 

6.0 

0.567 

0.201 

5 

15 

275.0 

2.59 

6.0 

0.6 

0.244 

6 

27 

400.0 

2.55 

6.0 

0.9 

0.1455 

7 

29 

300.0 

2.70 

6.0 

1.578 

0.1355 

8 

39 

51.75 

1.035 

6.0 

2.705 

1.082 

9 

44 

220.0 

5,43 

6.0 

0.9118 

0.227 

10 

48 

500.0 

7.808 

6.0 

0.74 

0.0595 

11 

64 

376.0 

7.14 

6.0 

0.532 

0.13 

12 

66 

127.8 

2.176 

6.0 

1.405 

0.486 

13 

68 

103.5 

1.473 

6.0 

2.28 

0.599 


Table 

G.2 

: Line Data 





Line 

No. 

From 

Bus 

To 

Bus 

Line 

R 

Impedance 

X 

— Y 

2 charge 

Turns 

Ratio 

1 

9 

8 

0.0 


0.0570 

0.0 

1.05 

2 

9 

7 

0.0320 

0.0780 

0.0090 

1.00 

3 

9 

5 

0.0660 

0.1600 

0.0047 

1.00 


Continued .... 
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Line 

Fxom 

To 

Bus 

Line 

R 

Impedance 

X 

1 Y 

^ chacge 

Turns 
Ratio . __ 

4 

9 

10 

0.0520 

0.1270 

0.0140 

1.00 

5 . 

10 

11 

0.0660 

0.1610 

0.0180 

1.00 

6 . 

7 

10 

0.0270 

0.0700 

0.0070 

1.00 

7 

12 

11 

0.0 

0.0530 

0.0 

0.95 

8 . 

11 

13 

0.0600 

0.1480 

0.0360 

1.00 

9 

14 

13 

0.0 

0.0800 

0.0 

1.00 

10 

13 

16 

0.0970 

0.2380 

0.270 

1.00 

11 

17 

15 

0,0 

0.0920 

0.0 

1.05 

12 . 

7 

6 

0.0 

0.2220 

0.0 

1.05 

13 

7 

4 

0.0 

0.0800 

0,0 

1,00 

14 . 

4 

3 

0,0 

0.0330 

0.0 

1.05 

15 , 

4 

5 

0.0 

0.1600 

0.0 

1.00 

16 

4 

12 

0 , Ol60 

0.0790 

0.0710 

1.00 

17 

12 

14 

0.0160 

0.0790 

0.0710 

1.00 

18 . 

17 

16 

0.0 

0.0800 

0.0 

0.95 

19 

2 

4 

0.0 

0.0620 

0.0 

1.00 

20 . 

4 

26 

0.0190 

O . D950 

0.1930 

1.00 

21 . 

2 

1 

O.o 

O.0340 

0.0 

1.05 

22 . 

31 

26 

0.0340 

0.1670 

0.1500 

1,00 

23 . 

26 

25 

0.0 

0.0800 

0.0 

0.95 

24 

25 

23 

0.2040 

0.5200 

0.0130 

1.00 

25 .. 

22 

23 

0.0 

0.0800 

0.0 

0.95 

26 > 

24 

22 

0.0 

0.0840 

0.0 

0.95 

27 . 

22 

17 

0.0480 

0.2500 

0.0505 

1.00 

28 . 

2 

24 

0.0100 

0.0200 

0.3353 

1.00 

29 . 

23 

21 

0.0366 

0.1412 

0.0140 

1.00 

30 . 

21 

20 

0.0720 

0.1860 

0.0050 

1.00 

31 

20 

. 19 

0.1460 

0.3740 

0.0100 

1.00 

32 , 

19 

18 

0 . 059C 

1 0.1500 

0.0040 

1.00 

33 

18 

16 

0.0300 0.0755 

0.0080 

1.00 

34 - 

28 

27 

0.0 

0.0810 

0.0 

1.05 

35 . 

30 

29 

0.0 

0.0610 

0.0 

1.05 


Continued 



296 


Line 

From 

Bus 

To 

.Bus 

Line Impedance 

. a _ X 

^ ^charge 

Turns 

Ratio 

36. 

32 

31 

0,0 0.0930 

0.0 

0.95 

37. 

31 

30 

0.0 0.0800 

0.0 

0.95 

38 

28 

32 

0.005 0.0510 

0.6706 

1.00 

39 

31 

33 

0.0130 0.0640 

0.0580 

1.00 

40 

31 

47 

0.0100 0.0790 

0.1770 

1.00 

41 

2 

32 

0.0158 0.1570 

0.5100 

1.00 

42 . 

33 

34 

0.0 0.0800 

0.0 

0.95 

43 

35 

33 

0.0 0.0840 

0.0 

0.95 

44. 

35 

24 

0.0062 0.0612 

0.2012 

1,00 

45.. 

34 

36 

0.0790 0,2010 

0.0220 

1.00 

46 

36 

37 

0.1690 0.4310 

0.0110 

1.00 

47 

37 

38 

0.0840 0.1880 

0.0210 

1.00 

48 

40 

39 

0.0 0.3800 

0.0 

1.05 

49. 

40 

38 

0.0890 0.2170 

0.0250 

1.00 

50 . 

38 

41 

0.1090 0.1960 

0.0220 

1.00 

51 

41 

51 

0.2350 0.6000 

0.0160 

1.00 

.52 

42 

41 

0.0 0.0530 

0.0 

0.95 

53. 

43 

42 

0.0 0.0840 

0.0 

0.95 

54 r 

47 

49 

0.0210 0.1030 

0.0920 

1,00 

55. 

49 

48 

0.0 0.0460 

0.0 

1.05 

56.. 

49 

50 

0.0170 0.0840 

0.0760 

1.00 

57 

49 

42 

O .037 O 0.1950 

0.0390 

1.00 

58- 

50 

51 

0.0 0.0530 

0.0 

0,95 

59 . 

52 

50 

0.0 0,0840 

0.0 

0.95 

60 V 

50 

55 

0.0290 0.1520 

0.0300 

1.00 

61 

50 

53 

0.0100 0.0520 

0.0390 

1,00 

62 

53 

54 

0.0 0.0800 

0.0 

0.95 

63 

51 

54 

0.0220 0.0540 

0.0060 

1.00 

64. 

55 

56 

0.0160 0.0850 

0.0170 

1.00 

65 

56 

57 

0.0 0,0800 

0.0 

1.00 

66. 

57 

59 

0.0280 0.0720 

0.0070 

1.00 

67 

59 

58 

0.0480 0.1240 

0.0120 

1.00 
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Line 

From To 
— Bus 

Line 
B.- 

Impedance 1 _ Turns 

X _2 a « +i- n 

66 

60 

59 

0.0 

0.0800 

0.0 

1,00 

69. 

53 

60 

0.0360 

0.1840 

0.0370 

1.00 

70. 

45 

44 

0.0 

0.1200 

0.0 

1.05 

71. 

45 

46 

0.0370 

0.0900 

0.0100 

1.00 

72 

46 

41 

0.0830 

0.1540 

0.0170 

1,00 

73 

46 

59 

0.1070 

0.1970 

0.0210 

1,00 

74 

60 

61 

6.0160 

0.0830 

0.0160 

1.00 

75 

61 

62 

0.0 

0.0800 

0.0 

0.95 

76. 

rs 

62 

0.0420 

0.1080 

0.0020 

1.00 

77. 

62 

63 

0.0350 

0.0890 

0.0090 

1.00 

78 

69 

68 

0.0 

0.2220 

0.0 

1.05 

79. 

69 

61 

0.0230 

0.1160 

0.1040 

1.00 

80. 

67 

66 

0.0 

0.1880 

0.0 

1.05 

81. 

65 

64 

0.0 

0.0630 

0.0 

1.05 

82. 

65 

56 

0.0280 

0.1440 

0.0290 

1.00 

83. 

65 

61 

0.0230 

0.1140 

0.0240 

1,00 

04 

65 

67 

0.0240 

0.0600 

0.0950 

1.00 

85 

67 

63 

0.0390 

0.0990 

0.0100 

1,00 

86. 

61 

42 

0.0230 

0.2293 

0.0695 

1.00 

87 

57 

67 

0.0550 

0.2910 

0.0070 

1.00 

88.. 

45 

70 

0.1840 

0.4680 

0.0120 

1.00 

89 . 

70 

38 

0,1650 

0.4220 

0.0110 

1.00 

90. 

33 

71 

0.0570 

0.2960 

0.0590 

1.00 

91.. 

71 

37 

0.0 

0.0800 

0.0 

0.95 

92 

45 

41 

0.1530 

0.3880 

0.0100 

1.00 

93. 

35 

43 

0,0131 

0.1306 

0.4293 

1.00 

94. 

32 

52 

0.0164 

0.1632 

0.5360 

1.00 



Tabid G.3 


Shunt Capacitor Data 


S,No* 

Bus No. 

Shunt Load 

Admittance 

1 

2 

0.0 

-0.4275 

2 

13 

0.0 

0.1500 

3 

20 

0.0 

0.0800 

4 

24 . 

0.0 

-0.2700 

5 

28 

0.0 

-0;.3375 

6 

31 

0.0 

0.2000 

7 

32 

0.0 

-0.8700 

8 

34 

0.0 

0.2250 

9 

35 

0.0 

-0.3220 

1C. 

36 

0.0 

0.100 

11 

37 

0.0 

0.3500 

12 

38 

0.0 

0.2000 

13 

41 . 

0.0 

0.200 




-0.2170 

14 

43 

W ♦ V 


15 

46 

0.0 

0.1000 

16 

47 

0.0 

0,3000 

17 

50 

0.0 

o.iooo 

18 

51 

0.0 

0.1750 

19 

52 

0.0 

-0.'2700 

20 

54 

0.0 

0.1500 

21 

57 

0.0 

0.1000 

22 

59 

0.0 

0.C75U 

23 

21 

0.0 

0.0500 






Table G.4 : Operating Data 


Bus 

Generation 


No * 


( mVaR ) 

r 

1 

630*667 

169.722 

0.0 

2 

0.0 

0.0 

0.0 

3 

506.0 

149.492 

0.0 

4 

0.0 

0.0 

0.0 

5 

0.0 

0.0 

0.0 

6 

98.0 

31.986 

0 . 0 , 

7 

0.0 

0.0 

12.8 

8 

297.0 

124.227 

0.0 

9 

0.0 

0.0 

185.0 

10 

0.0 

0.0 

80.0 

11 

0.0 

0.0 

155.0 

12 

0.0 

0.0 

0.0 

13 

0.0 

0,0 

100.0 

14 

0.0 

0.0 

0.0 

15 

184.08 

114.0 

0.0 

16 

0.0 

0.0 

73 . 0 - 

17 

0.0 

0.0 

36.0 

18 

0,0 

0.0 

16.0 

19 

0.0 

0.0 

32.0 

20 . 

0.0 

o:o 

27.0 

21 

0.0 

0.0 

32.0 

22 

0.0 

0.0 

0.0 

23 

0.0 

0.0 

75.0 

24 

0.0 

0.0 

0.0 

25 

0,0 

0.0 

133.0 

26 

0,0 

0.0 

0.0 

27 

304.0 

76.287 

0.0 


Power ■ ,„ Bu s . Vqi t . ag $— , 

( MVAR ) Mag .( pu ) Ang . ( deg ) 


0.0 

1.030 

0.0 

0.0 

1.057 

- 5.931 

0.0 

1.025 

- 1.406 

0.0 

1.054 

- 6.061 

0.0 

1.047 

- 5.064 

0.0 

1.025 

• 0,397 

8.3 

1.045 

- 5.723 

0.0 

1.025 

- 0.204 

130.0 

1.0433 

- 4.976 

50.0 

1.026 

- 7.413 

96.0 

1.023 

- 11.762 

0.0 

1.003 

- 10.057 

62.0 

1.007 

- 13.839 

0.0 

1.004 

- 11.984 

0.0 

1.005 

- 9.413 

45.5 

1.024 

- 16.682 

22.4 

1.005 

- 14.460 

9.9 

1.009 

- 17.881 

19.8 

0.981 

- 19.786 

16.8 

1.003 

- 21.838 

19.8 

1.016 

- 21.376 

0.0 

1.008 

- 17.625 

46.6 

1.036 

- 19.815 

0.0 

0.985 

- 16.193 

82.5 

1.039 

- 17.441 

0,0 

1.016 

- 14.382 

0.0 

1.025 

- 5.453 



Bus 

No . 


Generation 

■. Load Power 

_ Bus Voltage 

rMw)"" 

(MVrtRj 

( MW ) 

CMVAR ) Mag . ( pu ) Ang • ( deg ) 

28 

0,0 

0.0 

30.0 

20.0 

1.055 

- 12,342 

29 

261.0 

70.506 

0.0 

0.0 

1.025 

- 15.519 . 

30 

0.0 

0.0 

120.0 

74.5 

1.057 

- 19.947 

31 

0.0 

0,0 

160.0 

99.4 

1.014 

- 22.817 

32 

0.0 

0.0 

0.0 

0.0 

1.025 

- 115.936 

33 

0.0 

0.0 

0.0 

0.0 

0.998 

- 25.515 

34 

0.0 

0.0 

112.0 

69.5 

1.306 

- 29.639 

35 

0.0 

0.0 

0.0 

0.0 

0.966 

- 21.803 

36 

0.0 

0.0 

50.0 

32.0 

0.994 

- 34.030 

37 

0.0 

0.0 

147.0 

92.0 

0.967 

- 37.982 

38 

0.0 

0.0 

93.5 

88.0 

0.972 

- 37.638 

39 

25.0 

30.383 

0.0 

0.0 

1.025 

- 34.166 

40 

0.0 

0.0 

0.0 

0.0 

1.018 

- 36.913 

41 

0.0 

0.0 

255.0 

123.0 

liX57 

- 34.778 

42 

0.0 

0.0 

0.0 

0.0 

0.979 

- 31.484 

43 

0.0 

0.0 

0.0 

0.0 

0.946 

- 27.835 

44 

180.0 

55.037 

0.0 

0.0 

1,025 

- 24.448 

45 

0.0 

0,0 

0.0 

0.0 

1.048 

- 30.514 

46 

0.0 

OiO 

78.0 

38.6 

1.020 

- 32.678 

47 

0.0 

0.0 

234.0 

145.0 

0.989 

- 29.053 

48 

341.0 

256.0 

0.0 

0.0 

1.005 

- 25.781 

49 

0.0 

0.0 

295.0 

183.0 

0.997 , 

- 30.500 

50 

0.0 

0.0 

40.0 

24.6 

0.977 

- 32.056 

51 

0.0 

0.0 

227.0 

142.0 

1.005 

- 35.309 

52 

0.0 

0.0 

0.0 

0.0 

0.957 

- 26.552 

53 

0.0 

0,0 

0.0 

0.0 

0.972 

- 32.845 

54 

0.0 

0.0 

108.0 

68.0 

1.004 

— 35. 318 

55 

0.0 

0.0 

25.5 

48.0 

0.985 

- 31.313 

56 

0.0 

0.0 

0.0 

0.0 

1.013 

- 29.253 




57 

0,0 

0.0 

55.6 

35.6 

1.016 

- 29.967 

58 

0.0 

0,0 

42 <0 

27.0 

1,018 

- 30.261 

59 

0.0 

0.0 

57.0 

27.4 

1.013 

- 30,715 

60 

0,0 

0.0 

0,0 

0.0 

1.008 

- 29.974 

61 

0,0 

0.0 

0.0 

0.0 

1.019 

- 27.874 

62 

0.0 

0.0 

40.0 

27.0 

1.044 

- 28,940 

63 

0.0 

0,0 

33.2 

20.6 

1,042 

- 27.695 

64 

300.0 

72.882 

0.0 

0.0 

1.025 

- 19.401 

65 

0.0 

0.0 

0.0 

0.0 

1.057 

- 24.662 

66 

96.0 

25.636 

0.0 

0.0 

1.025 

- 20.610 

67 

0.0 

0.0 

14.0 

6.5 

1.055 

- 25.641 

68 

89.0 

26.689 

0.0 

0.0 

1.025 

- 19.829 

69 

0.0 

0.0 

0.0 

0.0 

1,050 

- 25.363 

70 

0.0 

0.0 

11.4 

7.0 

0,998 

- 34.790 

71 

0.0 

0.0 

0.0 

0.0 

0.927 

- 35.329 




o 


iO 


<n 






fig.g.1 single line diagram of 13-machine system 
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APPENDIX - H 


EL&nBsiTS OF [da/dx] AND [dg/dxi] IN EQ. (4,34) 


The elements of [da/dxj can be expressed as 


55^ 


Ej.V^ cos(©^ - 0^) 


i = 1,2 (m-2) 


(H-l) 


cos(9^ - 0^) 


1—1 ,2 , • • • • f (in“*2) 

^ — 1,2, ( in**l ) 

j 0 i 


(H-2) 




h cos(0|_-0|) 


1 — 1,2, ••••, ( m"*! ) 


(Ht3) 




m-2+i 


= 0 » 


1 = 1,2, ••••, 
j = 1,2, •••«, (m”*!) 

j ^ i 


(H.4) 


^92m- 


2m-2+i 


= 0 , 


i = 1,2, (n-m) 

j = 1,2, (m-1) 


(H.5) 



^^m-2+n+i 

Ml 


sin(©^ “ 


i = 1,2, .... , (m-l) 


(H.6) 


'in--2+n+i 


i = 1,2, ...., m 
3 = 1,2, ...., m-l 
3 ^ i 


(H.7) 


^^2m“»2-<-n+i 


'2n+m-l 




i=l,2, •••., (n-m) 
j = 1,2,'...., (m-l) 

i = 1,2, 


(H.8) 


(H.9) 


M. E. sin(©. - 0|) 


i = 1,2, ...-, 


(Hao) 


Mi E^. sin(G.| ~ 0.|) 


i = 1,2, («n-2) 

j =1,2, ...., ^ 
j ^ i 


(H.ll) 


'S'V. 

JL 


^ ~ + I B. . V. 

iii ^ 


sin 0ij 


df^A ) 


= ^ 


(H.12) 





*®J!=2±i 




1 ~ >•••} in 
j ~ ■IfSj •••») n 
3 ^ i 


(H.13) 


'2m-2+i 


jfx ®ij ^ij dV^ * 

i = 1»2» •.•«^ (n— in) 


(H.14) 


’2in-2±j 


^^m~24-n-fi 

w. 


V^ sin 0x j ^ ^ = 1»2» **••» (n—in) 

j ss Xf2> • m M • f n 

3 4 ^ 

2V^ Ej^ CQs(9j^ - 0 ^) 


(H,15) 


- 2BiiVi - 2 Vj cos 0ij 

01 ' 




+ * 

ay. 


i = 1,2, .... f m 


(H.16) 


*m-»24-n-fi 


^32m- 


2m-2+n+i 


- Vj^ cos 0x j * i — 1|2, ♦...,in 

j = 1,2, «...» n 

3 ^ ^ 


- ®lj''3 ’’iJ 


(H.17) 


df .(V.) 


i = 1,2, (n-m) (H.18) 



^92m-2+n+i 

m 


cos 0^^ » i = 1»2» ••••) (n— ni) 


j — •l-»2> ••••f n 

j i 


(Ha9) 


*2n4-in~l 


^ ^ 9 * f n 


(H.20) 


%Vi o°s(e^ - 0^) 

- M ; 


i = if2> 


(m-2) (H.21) 


W] 




m-2±i 


Mi E.V^. cos(9. - 0.) 


i = 1,2, «..., (m-2) 
j ~ 1,2, ••«•, n 


j 5^ i 


% Vj cos(0i - Qj) 


(H.22) 


+ Blj ViVj cos h} ’ 

i =: 1,2, ••••, m 


(H.23) 


^9m-2±i 


- Bij ViVj cos 0ij, 


^ ~ X^2f «•••! in 

j = X f2.f • • • • f n 

J 4 i 


(H.24) 
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^^2m-2+i 

= 


^92m-2+i 


^^m-2+n+i 

5PI 


^^m-2+n+i 

5^5 


^^2m-2+n+i _ 

Wl 


^^2in-.2+n+i 


n 

2 

j-1 


B, 


ij 


Vj 


1,2, 


• • » • I 


(n-m) (H.25) 


- %j ViVj cos 

i = 1,2, ..., (n-m) 

j — 1,2, ••••, n 
j i 

Ej^Vi sin(Gj^ - 0 ^) 


(H.26) 


n 


+ 2 B. . V.V. sin 0 . . 

Jssl -*• J "“"J 

X ss l|2f •#•*! ni 


(H.27) 


- hj Vj hi 


X lf 2 f ••••1^ in 
J = lf 2 f 0 • • 9 f n 

j ^ i 


(H.28) 


2 B^j V^Vj sin 0^y 

i = 1,2 (n-m) (H.29) 



V,iV- 
'^i J 


sin 0^j , 
i = 


1 , 2 , 


...,, (n-m) 


j = 1,2, ♦•••, n 

j 1 


(H.30) 



oua 


wi 


« 0 , 


1 5= lf2.t n (H«31) 


The elements of [dg/duj are as follows ; 


dg^ 

= 0, 

i ^ 1 f2. f ••••! 

(m-2) 

(H.32) 



cos (0j^ - ©j^) 





% 

i * k = m 


(H.33a) 


=s 0^ 

= 1^2) * » » » f 

k 0 i 

m 

(H.33b) 

^^2m-2+i 

«8k 

= 0, 

X =:'l;2f ••••) 

(n-m) 

(H.34) 

*®m-2+n+i 

%vk. 

sin(0j^ - ©j^) 






^ i = k = m 


(H.35a) 


= 0, 

i. — 1, f2. f ••••! 

k ^ i 

m 

(H.35b) 

^^2m-2+n+i 

= 0, 

1. =5 IfSy ••♦•f 

(n-m) 

(H.36) 

*92n+m-l 

55 ;; 

= 



(H.37) 
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.APPENDIX - I 

relationship between power injections in Lossy AND LOSSLSS 

NETWORKS [67] 


For an element ij (connected across nodes i and j in 
the network, the real and reactive power flows at node i, are 
given by 





cos 





sin 


(I.l) 



(vf - V^Vj cos 0ij) - g 






COS 


'ij 


( 1 . 2 ) 


From the above one can write 


^2 hg - rt? ''iJ 

wdiere a = 

Defining 

pij = '^ij 


(1.3) 


) (1.4) 


Q.' . 




(1.5) 

( 1 . 6 ) 


we have. 
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P ’ 


01 . 


Vk €C 

T77 « ■ IT? 

— ~ — 5 *’i1 ■'■ — ^ — 2 ®ii 


(1.7) 

( 1 . 8 ) 


Equations (1.5) and (1.6) describe a lossless network element. 
It is easy to see that the above equations still apply v\tien 
and Qjj^j are replaced by and and similarly for PJ^j 
and C^j. Thus, we have. 



1 

1 

-a 


h 

Si 

“ 1 . 2 

1 + a 

a 

1 


% 


(1.9) 


From Eq.(1.9), Eqs. (5.1) and (5.2) are easily derived. 



311 


APPENDIX - J 

H'JERGY STORED IN MACHINE REACTANCES AND TRANSMISSION ^JMES 

.Bjerqy Stored in Machine Reactances 

It was inentionad in Qiaptar 6 that ths energy stored 
in the machine reactances is given by 


m 

Z 

i=l 


^ vf - 2E^ V, cos(e, - 0^))]-^ 


- (cos(2S^ - 20j) - 1) 


V? 


(E^f + yf - 2Ei^ oos(e, . 0^)) 5^ . 5^] 

<?i qi 

(J.l) 

In general » it has been shown in [70] that the energy stored 
in a line element is equal to half the reactive power loss. 

It will be shown here that the energy term given in Eq.(J.l) 
is actually half the reactive power loss within the generator. 

The reactive power generation at the internal bus of 
X-th . machine is 


Qei = J - j = 


qi 


qi 


■dr 


F* i - E* i 
^di ^qi ^qi ""di 


(J.2) 
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The reactive power output at the terminal bus is 
Ji = I.n + 3 - J (J.3) 


Hence, the reactive power loss in the machine reactance is 


«ri = 

- ''di) - - \> 

(J.4) 

For the 

machine model considered, 



“ ^qi **■ ^di ^di 

(J.5) 

and 


(J.6) 

Thus, 

®ri ^ ^qi **’ 

(J.7) 

Also , 

^qi "" ^i 

(J.8) 


sin (©i - 0i) 

{J.9) 

Thus , 




^qi = e:» , + sin(Oi - 0^) 
cos(©i - 0^) - 

i^^ = -* -^:rr— » 


from Eqs.(6.12)and (6.l3) in Chapter 6 


(J.IO) 
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Hence* 


X* 

^qi ^qi 


V* 

^di ^di 


... r^di ■*■ 2 

■»* 

^ 

qi 

A ., ^°3(9, _- 0,) ■■ £ ^ 2 

^di ^ -» -* 


(Jai) 


Mi 


[Vj c os (9^ - 0^ - E« . ] 


(J.12) 


Mi 


‘<ri 


[e;,, -f V, sin(e^-0^)]^ ^ [V, cos(et-0,) - E'^f 


x“f 
""qi 


^di 


2 

Expanding the R.H.S. of Eq.(J.13) and replacing sin (0^ - 0^) 
by ^ [1 - cos 2(0^- 0^) and cos^(0j^ - 0^) by 
^ [cos 2(Q^ - 0^) + 1]* we get after rearranging the terms, 


^ri 


,2 

_ zsi 


E'r - 2E!. cos(0i - 0.) w >' n 

03. oi i i 1__ ^ j-^gg 2(0. -0.) + l] 

^Xdi 3. 1 


^di 


ELr + 2 E'.v. sin(0. -0.) r oc m 

+ -di ^ i 1 + [1 . cos 2{©^-0^)] 

(J.14) 


q3- 


V? 


Adding and subtracting the term — 7 — from the R.H.S* of Eq. 

""qi 

(J.I 4 ) and simplifying, we get. 
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[E>f + vf -2E'iViCoa(e.Ha.)J 


■V^ [cos(2>S,-20.)..l-] 2^ 

1 1 i 

V? 

■ [Eii W + 2E1, Vj sin(0,-0,) ^ - 4- ] 

qi qi 


(J.15) 


Eq*{J.15) gives the reactive power loss in the reactances of a 
single machine. The total reactive power loss in the reactances 
of machines for a system having m machines> is given by 


IB ni oo. 

Qxi = ^ ''i ■ ''l cos(ei-0i) )-^- 

- V2 (00.(20,-20,)-! 


+ (E^ + vf + 2Ei, Y, .ln(0,-0,)) ^ 


VI 


-] (J.16) 


qi 


This is tv/ice the energy term given in Eq.(J.l}- 


J - ^ ^nA -rgy fitored in the Tran smission Lines 

Energy stored in the transmission lines has been 

given by 

■ ^ ill jL 


(J.17) 



It will be shown tJiot the energy term given in Eq.(J.17) is 
half the reactive power loss in the transmission lines. 


Power injection at bus i of a lossless network can 
be calculated from 


Pi + 


Hence , 
Qi 


n j(0, - 04 ) 

B,j) V,V. e J 

n , . 

- j (cos 0ij + j sin 


= - hi ''j 


(Jel8) 


(J.19) 


In the above expressions - Im where [Yj is the 
admittance matrix of the transmission network (excluding 
machine reactances). 

The reactive power loss in the transmission network 
can be given by 


<3rl 

n 

* E CL = 

i«l ^ 

n 

- E 
i=l 

This 

is twice the 

energy 


n 


E V^Vj cos 0ij 


term given in (J.17) 


(J.20) 
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Thusf in view of Eqs.(J*16) and (J.20) the total loss 
in the stationary network (including generator stator windings) 
is twice the sum of the energy terms given in Eqs.(J.l) and 
(J.i7). The total energy stored is equal to, 


1 


5 


L 

2 

15=1 




(J.21) 


where* 




=5 reactive power generation at the internal bus of 
generator i 
=5 load at bus j 

= total number of elements in the network including 
machine reactances. 
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APPENDIX - K 

DATA FOR 4-MACHINE SYSTEM 


Number of 

machines 

* 

4 




Number of 

buses 

• 

• 

6 




Number of 

lines 

m 

* 

7 




Base MV A 


• 

100 




Table K,1 

* Machine 

constants 




Generator 

M 




^d 


No. 

in p.u 

2. 

sec /rad 


in pu 

in pu 

in pu 

in sec 

1 

0.5302 


0.004 

- 

•m 


2 

0.0159 


0.500 

0.950 

1*450 

8,0 

3 

0.0106 


0.400 

0*750 

1.150 

6.0 

4 

0.0080 


1.0 

- 

- 

- 
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Table K.2 s Line data 


Line No. 

From bus 

To bus Line reactance 

in pu 

1 

1 

5 


0,15 

2 

1 

4 


0.20 

3 

4 

2 

- 

0.50 

4 

2 

3 


0.80 

5 

3 

6 


0.30 

6 

6 

5 


0.40 

7 

4 

6 


0.50 

Table K,3 : 

Operating data 



Bus Voltaae 

Bus angle 

Generation 

Loads 

No, Magnitude (in deg.^ Active 

Reactive 

Active Reactive 

(in pu) 

power 

power 

power power 


in pu 

in pu 

in pu in pu 


1.0 

0.0 

0.30 

0.0875 


— 

1.0671 

5.760 

0.30 

0.20 

- 

- 

1.0263 

2.0263 

0.20 

0.10 

- 

- 

1.0048 

0.0938 

0.10 

0.05 

0.20 

0.10 

0.9843 

-2.6920 

- 

- 

0.30 

0.10 

0.9861 

-2.8850 

— 

— 

0.40 

0.15 


* 


6 
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